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PREFACE 

In the late 1960's the highly successful Mariner missions opened an 
exciting era of planetary exploration. Our knowledge of the near planets 
soon be further enhanced by the Mariner and Viking missions scheduled 
for the early 1970's. In the later half of this decade our exploration of 
the solar system will be extended to the outer planets - Jupiter, Saturn, 
Uranus , Neptune and Pluto. The alignment of the outer planets in the period 
from 1976 to 1980 offers several favorable opportunities for successive fly- 
bys of various combinations of the planets. One of the most interesting, 
and ambitious, is an opportunity in 1977 for launching a spacecraft on a 
mission which includes successive flybys of Jupiter, Saturn, Uranus and 
Neptune - the Grand Tour. 

Because of the extreme sensitivity of the Grand Tour mission to er- 
rors in the encounter trajectory at the intermediate planets, the success of 
the mission will depend to a large extent on the effectiveness of the navi- 
gation and guidance process during planetary approach. In this study, the 
theory of linear estimation was applied to the problem of defining the ap- 
proach navigation accuracy and guidance requirements for the Grand Tour 
mission. 
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ABSTRACT 


The navigation and guidance process for the Jupiter, Saturn and 
Uranus planetary encounter phases of the 1977 Grand Tour Interior mission 
was simulated. The objectives of the investigation were to (1) define ref- 
erence approach navigation accuracies, (2) evaluate the relative information 
content of the various observation types, (3) survey and define reference 
encounter guidance requirements, (4) determine sensitivities to assumed sim- 
ulation model parameters, and (5) assess the adequacy of the linear estima- 
tion theory. 

A linear sequential estimator is used to provide an estimate of the 
augmented state vector which consists of the six state variables of position 
and velocity plus the three components of a planet position bias. The dis- 
tribution of random observation errors is assumed to be Gaussian with zero 
mean and known covariance matrix. Navigation accuracies are determined from 
an analysis of the error covariance matrix associated with the state esti- 
mate. Evaluation of the errors in the estimates for the state and target 
vectors are accomplished by direct comparison with the simulated data. 

The guidance process is simulated using a non-spherical model of the 
execution errors. The uncertainty in the velocity correction is assumed to 
originate from: (1) uncertainty in the magnitude of AV due to a random error 
in the total impulse applied by the control rocket, and (2) uncertainty in 
the orientation of the correction vector due to random errors in the direc- 
tion of the thrust vector. A simplified guidance algorithm is developed and 
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used to compute the velocity correction required at a specified time to null 
the deviation between the estimated target vector and the nominal target 
vect or. 

Computational algorithms which simulate the navigation and guidance 
process were derived from the theory and implemented into two research-ori- 
ented computer programs. Programs STEP V and STEP VI were written in FORTRAN 
IV for the CDC 6600 system at the University of Texas at Austin. Program 
STEP VI utilizes the interactive graphics capability of the CDC 252 display 
console to provide a powerful computational tool for the rapid analysis of 
interplanetary navigation and guidance problems. 
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CHAPTER 1 


INTRODUCTION 


During the last half of this decade the United States will embark on 
one of the most significant scientific investigations of all time — the explor- 
ation of the outer planets. Unmanned spacecraft will travel billions of miles 
to the outer edge of our solar system and send back data which will help scien- 
tists understand the origin and evolution of the universe . The key to the in- 
vestigation will be the planet Jupiter, regarded by many as the primary object 
of interest in the solar system. Because of Jupiter 1 s powerful trajectory 
shaping capability , missions to the outer planets will include a swingby of 
Jupiter as a means of reducing launch energy requirements and associated trip 
times . 

The alignment of the outer planets in the period from 1976 to 1980 of- 
fers favorable opportunities for several multi-planet flyby missions , ^ 3 ^ ^ some 
of these are: 


1976- 80 Jupiter-Saturn , 

1978-80 Jupiter-Pluto , 

1978-80 Jupiter-Uranus-Neptune , 

1977- 78 Jupiter-Saturn-Pluto , and 

1976-80 Jupiter-Saturn-Uranus-Neptune (Grand Tour). 

Of particular interest is the four-planet flyby mission commonly referred to 
as the Grand Tour. The Grand Tour involves successive flybys of the planets 
Jupiter, Saturn, Uranus, and Neptune. 
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1.1 Mission Description 

The 1977 Grand Tour Interior Mission departs Earth in September of 1977 
and arrives at Neptune nine years later. A profile of the heliocentric trajec- 
tory for the nominal mission is shown in Figure 1. Jupiter encounter at an al- 
titude of 194,500 km occurs 550 days after Earth departure. The swingby of 
Jupiter turns the velocity vector and increases the heliocentric energy so that 
the spacecraft encounters Saturn 576 days later. At Saturn the spacecraft 
passes inside the rings at an altitude of 9,180 km. The Saturn-Uranus leg takes 
1230 days and brings the spacecraft within 13,100 km of Uranus. After an addi- 
tional flight time of 998 days the spacecraft swings by Neptune and departs on a 
trajectory which escapes the solar system. The heliocentric trajectory lies 
very close to the ecliptic plane throughout the mission with a maximum heliocen- 
tric inclination of 2.45 degrees for the Jupiter-Saturn leg. Total heliocentric 
angle displacement for the mission is 326 degrees. A more detailed definition 
of the nominal mission is given in Appendix A. 

1.2 Problem Description 

From a navigation and guidance standpoint, the Grand Tour mission pre- 
sents several unusual problems which have not been encountered in previous in- 
terplanetary missions. Some of these are listed below: 

(1) Severe guidance requirements . Extreme trajectory error sensitivi- 
ties result in a requirement for guidance maneuvers during approach 
and during departure at each planet . ^ 5 ^ ^ 

(2) Long round-trip communication time. In addition to operational 

problems caused by the time delay (5.5 hours at Uranus), the ac- 

(7) 

curacy of the earth-based doppler data is reduced. 



NEPTUNE ENCOUNTER 
V 9/26/86 (9 YEARS) 



/ 5AU' 

/ / 

EARTH LAUNCH 

I 9/4/77 


URANUS ENCOUNTER 
T\l/2/84 (6.25 YEARS) 


I > SATURN ENCOUNTER 
Ls\ 8/20/80 (3 YEARS) 

JUPITER ENCOUNTER 
1/22/79 (1 .25 YEARS) 


Figure 1 Heliocentric Profile of the 1977 Grand Tour Mission 
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(3) uncertainties in the planetary ephemerides . Uncertainties 
on the order of 0.2 arc-sec are predicted for the late 1970 's^ 8 \ 
(At Uranus 0.2 arc-sec is equivalent to about 3,000 km.) 

(4) En counter constraints. The trajectory must either pass outside 
Saturn's rings (exterior mission) or through the. small region in- 
side the rings (interior mission). Exterior missions have longer 
total mission durations, lower launch energy requirements, and 
less critical guidance requirements than the interior missions. 

(5) Long mission durations. The long trip times (typically about 9 
years) and extreme trajectory sensitivity make precision numeri- 
cal simulations both difficult and expensive. A high degree of 

reliability is required for onboard navigation and guidance sys- 
terns. 

Because of the extreme sensitivity of multi-planet flyby missions to 
errors in the encounter trajectory at the intermediate planets, the success of 
the Grand Tour will depend to a large extent on the effectiveness of navigation 
and guidance procedures during planetary approach. The spacecraft will arrive 
at the sphere of influence of each planet with large uncertainties in the state 
of the spacecraft relative to the planet. Approach navigation must reduce this 
uncertainty to an acceptable tolerance before an effective guidance correction 
can be performed. Guidance corrections based on poor navigation information 
will result in an incorrect execution and wasted guidance propellant. If the 
execution errors remain uncorrected throughout approach they will be amplified 
by the encounter and will require large post-encounter corrections. 

Approach guidance correction requirements vary with execution time and 
increase sharply a few days prior to encounter. Navigation accuracy also in- 
creases during approach due to the improved planetocentric state information 
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provided by the observations. Since increased approach navigation accuracy is 
reflected in lower post-encounter corrections, the selection of the approach 
guidance execution time for minimum guidance propellant is an important opti- 
mization problem. 

1.3 Technical Approach 


In this study the navigation and guidance process for the Jupiter, 
Saturn and Uranus planetary encounter phases of the 1977 Grand Tour Interior 
mission has been simulated. The objectives of the investigation were to (1) 
define reference approach navigation accuracies, (2) evaluate the relative in- 
formation content of the various observation types, (3) survey and define ref- 
erence encounter guidance requirements, (4) determine sensitivities to assumed 
simulation model parameters, and (5) assess the adequacy of the linear estima- 
tion theory. 

A linear sequential estimator is used to provide an estimate of the 
augmented state vector which consists of six state variables of position and 
velocity plus the three components of a planet position bias. The distribu- 
tion of the random observation errors is assumed to be Gaussian with zero mean 
and known covariance matrix. Reference navigation accuracies are determined 
from analysis of the error covariance matrix associated with the state esti- 
mate. Evaluation of the accuracy of the estimated state and target vector are 
accomplished by direct comparison with simulated data. 

The guidance process is simulated using a non-spherical model of the 
execution errors. The uncertainty in the velocity correction is assumed to 
originate from: (1) uncertainty in the magnitude of the correction vector due 

to a random error in the total impulse applied by the control rocket, and (2) 
uncertainty in the orientation of the correction vector due to random errors 



in the direction of the thrust vector. A simplified guidance algorithm is 
developed and used to compute the velocity correction required to null the 
deviation between the estimated target vector and the nominal target vector. 

In Chapter 2 the computational algorithms and error models used to sim- 
ulate the navigation and guidance process are derived from the theory. Sup- 
porting derivations are found in Appendices B, C, and D. The computational 
algorithms were incorporated into two research-oriented computer programs. 

The programs, one an interactive graphics program, are described in Chapter 3. 

Chapter 4 contains a comprehensive analysis of the navigation process 
during approach to Jupiter, Saturn, Uranus and Neptune. Reference navigation 
accuracies are defined for several combinations of Earth-based radar and on- 
board angle observations. Appendices E and F supplement Chapter 4. 

An analysis of approach guidance is contained in Chapter 5. Guidance 
execution times are surveyed to determine the variation of guidance require- 
ments and target accuracy during approach to Jupiter, Saturn, Uranus and 
Neptune. Sensitivities are defined for the effect of observation accuracy, 
guidance execution accuracy, and initial state errors on the guidance process. 
Finally , the results of a Monte Carlo simulation of the approach navigation 
and guidance process at Saturn encounter are presented. Appendices E, F and 
G supplement Chapter 5. 

Chapter 6 summarizes the pertinent results and conclusions of the 
study and outlines recommended extensions to the work. 



CHAPTER 2 


THEORETICAL DEVELOPMENT 


This chapter contains the theoretical development of the computational 
algorithms used to simulate planetary approach navigation and guidance. The 
chapter is divided into five sections as follows: 2.1 State Equations, 2.2 Ob 

servation-State Relationships, 2.3 Error Simulation Models, 2.4 Linear Estimation 
Theory, and 2.5 Guidance Theory. 

2.1 State Equations 

In the analysis presented here, the spacecraft motion relative to the 
target planet is approximated by a three-dimensional elliptic restricted three- 
body model. For this model the target planet is assumed to move about the sun 
in an unpreturbed elliptic orbit while the third body (the spacecraft) moves 
under the influence of both primaries. The state (position and velocity) of 
the spacecraft relative to the target planet is obtained by numerical integra- 
tion of the total acceleration resulting from the attraction of the sun and 
the planet. 

Equation of Motion 

The motion of a spacecraft which is perturbed by the sun as it ap- 
proaches a target planet, is governed by the vector differential equation 

r = V(U + U, ) (2.1) 

o 1 


7 
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where r is the position vector of the spacecraft relative to the planet, 

U q = ~ is the two-body force function, y is the gravitational parameter of the 
planet, r is the magnitude of r, and is the third body disturbing function. 
The disturbing function, which expresses the perturbing effect of the sun on 
the spacecraft, is given by 



( 2 . 2 ) 


where y is the gravitational parameter of the sun, R is the magnitude of the 

heliocentric position vector of the target planet, R , and R is the magnitude 

of the heliocentric position vector of the spacecraft, R. Recognizing that 
1 r 

V(~.) = -(— ) 5 the equations of motion can be written as 
r 



The vector r, R and R^ are defined in Fig. 2. 

Eq. (2.3) can be reduced to a system of first order differential equa- 
tions by the following transformation 


r - v 



The heliocentric position of the planet, i.e., R^_ in Eq. (2.4), is obtained 
from a classical solution of the two-body problem. 


The State Vector 

The two first order, nonlinear, vector differential equations (Eqs. 
(2.4)) can be represented by 





^ A 

X = F(X,t ) 


where the state vector X is defined as 
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(2.5) 


% 

X = 


( 2 . 6 ) 


and the 6-vector F is the corresponding right hand side of Eqs. (2.4). 

The Augmented State Vector 

It is sometimes necessary or desirable to estimate certain parameters 
of the dynamical and/or error models whose values are not well established. 

For the problem treated here, the uncertainty in the position of the outer 
planets was considered large enough to justify its inclusion as a vector to be 
estimated. This was accomplished by augmenting the state vector with the 

planetary bias vector b. The augmented state vector, a 9-vector, is defined 
as 

r 

X = v . (2.7) 

_b_ 

Since the position vector of the target planet changes very slowly during the 
relatively short encounter period, the error in the position vector was as- 
sumed to be a constant bias. The augmented state vector is therefore governed 
by the vector differential equation 


X = 


F$,t) 


* J 


( 2 . 8 ) 


Analytical expressions for the components of Eq. (2.8) are given in rectangu- 
lar cartesian coordinates in Appendix B. 
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Linearization of the State Equations: The State Transition Matrix 

As discussed in the previous section the n-dimensional state vector, 
X, can be defined as the solution to the first order, nonlinear, vector dif- 
ferential equation 

X(t) = F(X,t) ; X(t ) = X . (2.9) 

o o v / 

Expanding Eq. (2.9) in a Taylor series about a nominal solution (designated 
with an asterisk) yields 


X = X* + I3F/3XJ* (X - X*) + ... 


( 2 . 10 ) 


Neglecting terms of order higher than the first Eq. (2.9) can be written 


x = Ax • x ( t ) — x 

o o 


( 2 . 11 ) 


where x = X - X* and A = [3F/3X]*. 

By the theory of linear differential equations ^ 9 ^ the solution of Eq. (2.11) 
can be expressed in the form 


x(t) = $(t,t )x 
o o 


( 2 . 12 ) 


where *(t,t Q ) is an nxn transition matrix relating state deviations at time 
to to state deviations at time t. Differentiating Eq. (2.12) yields 


x(t) = $(t,t )x . 

o o 

Subsequent substitutions result in the following development 

*Ct,t o )x Q = ACt) x(t) 

3>Ct,t )x = ACt) $Ct,t )x 
o o o o 

0 = [A(t) $(.t,t Q ) - $(t ,t Q )J x o . 


(2.13) 


(2.14) 
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For a nontrivial solution x will be non-zero. Therefore Eq . (2.4) can be 

o 

satisfied only if the matrix differential equation governing the state tran- 
sition matrix satisfies the relation 

i(t,t ) = A(t ) <Kt,t ); $(t ,t ) = I . (2.15) 

o o o o 

The state transition matrix has the following properties: 

1) $(t,t) = I, the nxn identity matrix 

2) $(t. ,t . ) = $(t. ,t. ) for any t., t. and t, 

k j k i i ] i j k 

“1 

3) $(t.,t.) = 0 (t.,t.) for any t. and t. 

1 J 3 i i 3 

For this study the transition matrix was obtained by numerical integration of 
the system of linear differential equations given by Eq. (2.15). Alternately, 
the nonlinear differential equation governing the state (Eq. 2.9) can be in- 
tegrated n times from t^ to t to generate the required matrix of numerical 
partial derivatives. The transition matrix is then formed by 

* (t >V = wrr • (2 - 16) 

o 

Analytical expressions for the elements of the 9x9 A matrix are given in 
Appendix B. 

2 . 2 Observation-State Relationships 

At any given time the observation vector Y is related to the state 
vector X by the nonlinear equation 


YCt) = G(X,t) + e 


(2.17) 
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where e is the random error in the observation. Expanding Eq. (2.17) in a 
Taylor series about the nominal state and neglecting terms of order higher 
than the first yields 

Y(t) = Y*Ct) + [3G/3XJ* (X(t) - X* ( t ) ) + e . (2.18) 

When this linearization is valid, the deviation in the observation can be 
expressed as 

y (t ) = H(t ) x(t ) + e . (2.19) 

where y=Y-Y*,x=X- X* and H(t) = J3G/3XJ*. 

Recalling from the previous section that 

x(tj.) = *(t iS t k ) x(t k ) (2.20) 

the observation deviation at time t^ can be expressed in terms of the state 
deviation at time t^ by 


y(' t i ) = H(t i ) 3>(t i ,t k ) x(t k ) + e i . (2.21) 

It is convenient to write Eq. (2.21) as 

y i = H i \ + £ (2.22) 

&nd this convention will be used in the following sections . Note that when 
i = k, then ^ . Analytical expressions for the elements of the ft matrix 

in rectangular cartesian coordinates are given in Appendix C. 

Observation Equations 


Four observation types were considered in this study: 
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1 ) Range (p ) 

2) Range-Rate (p) 

3) Sun-Planet Angle 

4) Star-Planet Angle 



Earth-based (Fig. 3) 


onboard (Fig* 4) 


When all four observations occur simultaneously, Y is a 4-vector defined as 
follows 

Y T = Ip, p, a, SJ. (2.23) 


The Range, the linear distance from the tracking station to the spacecraft, 
is defined by 

p = (p-p) 3 * (2.24) 


where p = R - R ,R is the heliocentric distance of a topocentric tracking 
s s 

station, and R is the heliocentric distance to the spacecraft. 

The Range-Rate , the time rate of change of range, is defined by 

p = (p*p)/p (2.25) 

where p = R - R . 

s 

The Sun-Planet Angle , the smaller angle between the spacecraft-planet line 
and the spacecraft-sun line, is defined by 


a 


-1 

cos 


r-R 
r R 


(2.26) 


or 


a 


-1 

cos 


r*R^ + r 
?~R 



(2.27) 


where 0 < a < tt and where r is the planetocentric position vector of the 
spacecraft . 
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The Star-Planet Angle , the smaller angle between the probe-planet line and the 
line from the probe to a navigation star, is defined by 

6 = cos ^[-(r*s)/rj (2.28) 

where s is a unit vector in the direction of the navigation star. 

Motion of the Tracking Stations 

Range and range-rate measurements were assumed to be taken from a 

tracking station whose heliocentric position vector R is given by 

s 


R = R + r (2.29) 

s e s 


where R is the heliocentric position of the earth, and r is the geocentric 

6 S 

position of the tracking station. Given the right ascension, a , declina- 

s 

tion, 6 , and geocentric distance, r , of the tracking station at time t, the 
s s 

heliocentric position of the station is defined by 


R 

s 


R + [Ejr 
e s 


cos 6 cos a (t ) 
s s 

cos 6 sin a (t) 
s s 

sin 6 

s 


(2.30) 


where [E] is a 3x3 rotation matrix which transforms the right-ascension de- 
clination system to the heliocentric system. If the heliocentric system is 
the ecliptic system, [E] is defined as 


[E] 


1 0 
0 cos i 

e 

0 -sin i 

e 


0 

sin i 

£ 

cos i 

£ -J 


(2.31) 


The argument i is the obliquity of the ecliptic. Tracking station geometry 
is illustrated in Fig. 5. 
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Figure 5 Tracking Station Geometry 
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With the assumption that the earth's rotation is uniform, the right 
ascension of the tracking station can be expressed as 


a <t) = a (t Q ) + a g (t - t Q ) (2.32) 

where a (t ) is defined at some epoch time t . Differentiation of Eq. (2.30) 
so ° 

with respect to time yields 


R 

s 


R + {E]r a 
e s s 


-cos 6 sin a (t) 
s s 

cos <$ cos a (t) 
s s 

0 


(2.33) 


2.3 Error Simulation Models 

Errors in the components of the observation, initial state and guid- 
ance velocity correction vectors were modeled as normally distributed scalar 
random errors with zero mean and known variance. The errors were simulated 
by sampling at random from a standard normal distribution (zero mean and unit 
variance) and then scaling the error by the square root of the given variance. 

The density function for the standard normal error z is given by 


-h / 2 


p(z) = (2tt) 2 exp(-z /2) . 


(2.34) 


The probability that a sampled value of the standard distribution will be less 
than a given value z is given by the cumulative distribution function^ 10 ^ 


F = ( 2tt )" 


exp(-z /2)dz 


(2.35) 


( 11 ) 

The inverse of Eq. (2.35) can be approximated by the curve fit equation 


z = T - 


c + c.r + c 0 r 2 

O J. Z 

l + d 1 r + d 2 r2 + d 3 r^ 


(2.36) 
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T = [ln(F _2 )] 2 (2.37) 

where the coefficients and d i have the following values: 

c = 2.515517 d = 1.432788 

0 1 

c = 0.802853 d. = 0.189269 

1 l 

c = 0.010328 d = 0.001308 . 

3 J 

Sampling of the standard normal distribution was accomplished by 
first sampling at random from a uniform distribution to obtain a value for 
F (0 £ F <; l)and then computing the standard normal error from Eq. (2.36). 
The simulated error is then computed as the product of the standard error 
and the square root of the given variance. 

Simulation of the Observations. 

The observational data were simulated by adding random errors to the 
observation value computed from the simulated value of the true state, i.e.. 


Y 


+ 

true 


z . 


(2.38) 


The true state is obtained by numerically integrating Eq. (2.4) using the 
simulated initial state vector. If the components of the error vector are 
assumed to be uncorrelated, the observation covariance matrix can be written 


R 



p J 


(2.39) 


0 


0 


0 



20 


The errors in the components can then be generated independently by the method 
of the preceding paragraph. 

Simulation of the Initial State Errors 

For the purposes of mission simulations, an initial state error is 
specified in one of two ways: 1) input an initial state error, or 2) gener- 

ate a random state error. When the second method is used, the components 
of the error vector are assumed to be uncorrelated with each other and are 
generated by the method described in the previous paragraphs. 


Simulation of the Guidance Correction Errors 

The uncertainty in the commanded velocity correction AV was assumed 
to originate from two independent sources: l) uncertainty in the magnitude 

A 

of AV due to a random error in the total impulse applied by the control rock- 
et, and 2) uncertainty in the orientation of the correction vector due to 

random errors in the direction of the thrust vector. The error model adopted 

( 12 ) 

for this study and derived below was first proposed by Battin 

Consider a coordinate system with the z-axis oriented along the com- 
manded correction vector AV as shown in Fig. 6. In the reference coordinate 
system the correction vector is given by 


AV 


AV D 


(2.40) 


where D is the transformation matrix which relates the new coordinate system 
to the reference coordinate system. 

With the assumption that the error in the magnitude of the correction 


is a proportional error, the magnitude of the actual correction can be written 

| AV I = (1 + k)|AV| (2.41) 





22 


where k is a small random variable. The actual correction vector is then 


given by 


sin y cos 0 


AV= (1 + k) AV D sin y sin 0 


cos y 


where Y is a small positive random variable and 0 is a random variable which 
is uniformly distributed in the interval -ir to tt . The error parameters k, 

Y 5 end 0 are assumed to be statistically independent. The error in the ve- 
locity correction vector can be expressed as 


n = AV - AV . 


With the small angle assumption for y, Eq. (2.43) can be written 


n = | AV | D <(1 + k) y sin 0 


cos 0 


+ K 0 


Now consider the covariance N of the error vector n, i.e., 


N = cov [n ] = E{nn ] . 


Substitution of Eq. (2.44) into Eq. (2.45) yields 


N = E [ k 2 .] E [AV AV ] t E[y 2 J(1 + x 2 ) D 0 1 0 D T . (2.46) 


10 0 


0 0 0 


Neglecting the term involving k 2 and simplifying, Eq. (2.46) becomes 


N = c 2 V + _L 
x 2 


V = E{AV AV ] 


{(Trace V)I - Vj 


(2.47) 
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where k and y are assumed to be distributed with zero mean and known variances 
° K Z = E[k 2 J> and a^ 2 = E[y 2 ]. 

2,4 Linear Estimation Theory 

When the dynamical system and the observational relationships can be 
satisfactorily approximated by linear equations, Eqs. (2.12) and (2.19) re- 
spectively , the well developed theory of linear estimation can be employed to 
determine the "best” estimate of the state. The statistical principle of max- 
imum likelihood is used here to derive a sequential estimation algorithm for 
interplanetary navigation. The procedure followed is described in greater de- 
tail in Reference 13. Under the assumptions of this study the resulting al- 
gorithm is the same as that obtained using the principal of minimum variance. 

The estimation problem considered here can be stated as follows. Find 
the maximum likelihood estimate of the epoch state given the sequence of 

observation vectors y^, Y 2 »***»y which contain random errors. Because the 
^ n Etial state (injection condition) x^ is unknown and because the mathematical 
model used to describe the dynamical system is inexact, the true epoch state 
x^ is never known. 

Maximum Likelihood Estimation 

Consider a set of m indpendent random samples y , y , . . . ,y drawn 

12 J m 

from a population characterized by the joint probability density function 

p(y 1 ,y 2 ,...,y ni ;x) = p^y^x) p 2 (y 2 ,x) ... P m (y m ,x) 

where x is a parameter of the distribution. The likelihood function is de- 
fined as 

m 

L = * P-. Cy. ,x) . (2.48) 

i=l 1 1 
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The value of x which maximizes Eq. (2,48) is defined as the best estimate x 
(in the maximum likelihood sense) and is called the maximum likelihood esti- 
mate of x. 

The Maximum Likelihood Principle may be stated as follows : select x 

to maximize the probability of the occurrence of the sample actually observed 
That is to say, of all possible samples the one that was actually observed is 
the one most likely, and therefore x should be chosen to maximize the probabil 
ity of its occurrence. 

Statistical Assumptions 

Consider the m linear observation-state relationships 

y. = H.x n + e. ; i = l,2,...,m (2.49) 

J i l k x 

where , the observation error vectors, are assumed to be normally distri- 
buted about a zero mean and to have a known covariance R, i.e., 

E[e.J = 0 , Ele^d] = R5„ ; i,j = 1,2,. . . ,m (2.50) 

where 6.. is the kronecker delta, i.e. , 6. . = 0, if i / i and 6.. = 1, if 
13 13 13 

i = j. Hence, it is assumed that the observation errors are not correlated 
in time , i.e., 

E ^ G i e j T ^ = ’ i / j * (2.51) 

. % 

The injection condition, x , and the error in its estimate, x = 

o o 

x q - x^, are assumed to have the following apriori statistics: 

E[x ] = E[x 1 = 0 
o o 

0/ % T T 

E [x x J = E{x x ] = P . 
o o o o o 


(2.52) 
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It is further assumed that the injection condition, x q , and the ob- 
servation errors, are not correlated, i.e.. 


EUq ] = 0 , for all k. 


(2.53) 


The Likelihood Function 

The best estimate of the actual state deviation vector will be the 
vector x^ which maximizes the likelihood function 


L = p( >v y 2 v V = p(y i y ,JV 


(2.54) 


where p(y,. ,y 0 > • • • ,y ,x, ) is the joint density function of observations and 
1 A m K 

state deviation x^ and where p(x^) is the density function of the random 
state vector x^. If no apriori statistics on x^ are considered then p(x^) 
is constant and the problem reduces to maximizing the conditional density 
function p (y 1 ,y 2 , • • . ,y m | \ ) • 

Given the multivariate Gaussian density function for the observation 
error vector e . 


-p/2 ■ -- 1 - 

Ir i' 


p(e.) = ( 2 tt ) t " | R. | ^ exp[-^ eT Rjt ep 


(2.55) 


and introducing the transformation 


e i = y A - H i*k 


H.: 


(2.56) 


the following conditional density function is obtained 


p ( y i = ( 2lT ) 


-p/2 


R 


T„-l> 


i ' 2 expt-35(y i - H^) R. (y £ - H.^)]. (2.57) 


With the assumption that the observation vectors are uncorrelated in 
time, the conditional density function can be written 
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m 

P(y l‘ y 2 V = P<y i V 


pm 

2 m 


(2.58) 


m 


= ( 2tt ) z tt | R | ^ expI-Jg y , (y. - H x ) T R 1 (y - H.x )] 
i=l 1 i=l 1 1 K 1 l l k 


Assuming that the true state is normally distributed about a mean 
conditioned on apriori information only, the density function for x^_ can be 
written 


P^) = 2tt n/2 |P k | * exp[-^(x k - x k ) T P k 1 (x k - x^ ) 3 (2.59) 

where x. = $(t, ,t )x is the mean at time t. , x = x is the apriori mean, 
k k ? o o k o o r 

T 

and P = $(t, ,t )P $(t ,t ) is the apriori state covariance matrix mapped 
K K O O K O 

to time t 

k 

From Eqs. (2.53) and (2.54) the likelihood function Eq. (2.49) can 
be written 


L = <2„r (mp+n)/2 




R. 2 (y . -H. x. 
i •'l i k 


) 


+ <x k - V T V 1( \ - *k ,I} • 


(2.60) 


Maximizing the Likelihood Function 

Since log L is a monotonic function and reaches its maximum value 
when L is a maximum, maximizing log L is equivalent to maximizing L. For 
convenience define 


J = log L = J - \ J 2 


(2.61) 


where 
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h = log [( 2 .)- ( " P +n)/2 | R|-"'' 2 | P k |^, ( 2 . 62 ) 

IT) T A 

J 2 = ( ^i - H i x k ) R i (yi - H i\) 

( 2 . 63 ) 

+ ( \ - \> T V 1 <\ - \> • 

For J to be a maximum with respect to the choice of ^ requires that 
the term be a minimum. For an extremal solution the first variation of 
J 2 must vanish, i.e., 

^ T _ 1 

6J = 0 = Z [-(y. - H,3t r R. H. fix. 

z i=l i x k ill 

- (H. S ^) 1 rT 1 (yj - H.^)] ( 2 . 64 ) 

t<i k- i k )I V 1 ' _1 

where x^ is that value of x^ which maximizes L. 

Since the first two terms on the right hand side are the scalar trans- 
pose of each other, they are equal. Likewise the second two terms are equal. 
Since <5x^ is arbitrary, it is necessary that 
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To insure that the solution obtained is a minimum the second varia- 
tion of with respect to x^ must be positive definite, i.e., 


rp in rn y| _ i 

° 2J 2 - S \ It, h H i + V 1Sx k ’ 0 ' (2 ' 67) 

Therefore, for Eq. (2.66) to maximize the likelihood function, Eq. (2.67) 
must be satisfied. 


Batch Processing 

Equation (2.66) can be used to compute the best estimate of the epoch 

state, x, , given the observation vectors y y and apriori statistics x 
k 1 m o 

and P^. The procedure is called ”batch processing” if all the observations 
are processed in a single batch. Several points regarding Eq. (2.66) should 
be noted before proceeding: 

1) The solution requires the processing of all observational data. 

2) The solution requires the inversion of an nxn matrix (first 
bracketed matrix in Eq. (2.66)). 

3) If no apriori data is considered, at least n independent scalar 
measurements must be processed to insure that the matrix to be 
inverted is full rank. 

4) If apriori information is available, a single scalar observation 
can be processed to provide a new estimate for all n components 
of the state vector. 

5) After processing all observational data to obtain a solution for 

the initial state x^, a new nominal trajectory may be generated 

and the data reprocessed to obtain an improved estimate for x . 

o 

Successive iterations will reduce errors caused by nonlinear ef- 


fects . 
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Sequential State Estimation 

For the special case when a single p-vector of observations is pro- 

th 

cessed (i.e., m = 1) and the i observation coincides with the epoch time 
t^ , Eq. (2.66) reduces to 

= V R_1 “k + V 1 * 1 V R_1 ^ + V 1 • <2 - 68) 

It can be shown in a straightforward manner that the first term on 

the left hand side of Eq. ( 2.64) is the covariance matrix of the estimate 

a, 

error x = x - x , i . e . , 


p k ° E[x k - - s k > ] 

-1 

T -1 _ _ 1 

= <H k r «k + p k > 


(2.69) 


Note that x^ and can be obtained by mapping and P ^ through the 


(9) 

transition matrix , i.e. , 


\ 3 ,I(t k’ t k-i )i k-i 


(2.70) 


P k $(t k 5t k-l )P k-l * ^k^k-l* 


(2.71) 


Thus Eq. (2.68) can be used to sequentially estimate the state after each ob- 
servation. Note, however, that the solution of Eq. (2.68) requires the inver- 

- -1 

+ ) . The order of the matrix inver- 

sion required to compute can be reduced from nxn to pxp by application of 
the Schur identity to Eq. (2.68). The result, derived in Appendix D, is 

P k = P k ' P k (H k P k H k T + R) 1 H k P k * (2.72) 

Eq. (2.72) offers two important computational advantages over Eq. (2.69). 
First, the matrix inversion required is of the order pxp rather than nxn. 


sion of the nxn matrix, (H^ R 
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In particular, if a single observation is made (p=l), Eq. (2.72) requires only 
a scalar division. Secondly, Eq. (2.72) does not require the inverse of ma- 
trix P. as does Eq. (2.69). 
k 

It will be convenient in the discussion which follows to define the 
nxp matrix 

x* = K \ T (H k H k T + R)_1 • (2 ' 73) 

Eq. (2.72) can then be expressed as 

p k - a \ v f k • <2 - ,4) 

Substituting Eq. (2.74) into Eq. (2.68) leads to 

\= U-h (H k T X ' 1 y k + V 1 *k> <2 ' 75> 

which can be rearranged to give 

\ = w - \ H k )5 k + (I - 1 Vh "k T R_1 y k • <2 - 76) 

The coefficient of y k in Eq. (2.76) is just as is shown in the following 
development: 

a - \ VW * 1 - 'W * 1 - f A I<H k h H k T + Rrl VA V ‘ 

- W 1 * 1 - ( Vk«k T + Rrl VA V ‘i 
= f k H k T<H k ? k H k T 7 R)_1 1 (H kVk T + R)R_1 - Vk"k T R ' 1] (2 - 77) 

■ f a t <va t R R) ' 1 
= 1 
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Thus Eq. (2.71) can be written 

x k = *k K k H k\ + K k y k 
which can be rearranged to give 

\ = \ + ‘ 


(2.78) 


(2.79) 


Sequential Estimation Algorithm 

The computational algorithm for the sequential estimation scheme de- 
rived here can be summarized as follows : 


1 ) 

2 ) 

3) 

4) 

5) 

6 ) 


Given the observation at the current time t^, and x^_ 1 and P^^ 
at the previous observation time t k _^ : 

Integrate X = F(X,t) from to t^ with initial condition X(t k _p 


= X*(t k _ 1 ) to obtain X*(t k ) 


Integrate $ = A$ from t k _^ to t k with initial condition 4> k _^ - I 


to obtain $(t k ,t k _p 


Compute Y * = G(X k *, t k ) and evaluate y k = Y k - Y k * . 
Evaluate H k using X*(t k ). 

Compute and P k from 


K = ‘'VWVi 


P k " 4 ’ <t k’ t k-l ,P k-l * 'Wd' 


7 ) Compute from 


■S: = U T (H k ? k C + R)_1 ' 


8) Compute x k and P k from 



\ = \ + “k&k - “k *k ] 

p k = t 1 - \ vv 


9) Store and . 


10) Repeat the process at the next observation time t, „ 

k+1 


2 . 5 Guidance Theory 

The purpose of the guidance system is to provide the necessary control 
of the velocity vector to achieve specified mission objectives. For this 
study the mission objectives are specified in terms of a 3-vector of target - 
ing parameters . The impulsive velocity change required to null the predicted 
deviation in the targeting vector is expressed as a function of the state de- 
viation. With the usual assumption of small state deviations, a guidance law 
is derived which expresses the commanded velocity correction AV as a linear 
function of the state estimate x. This linearization makes possible the de- 
velopment of a simple computational algorithm for planetary approach guidance. 

Targeting Vector 

In general , the targeting vector T is some nonlinear function of the 
state, i.e. , 

T = T(r,v). (2.80) 


For this study the targeting vector was defined as 
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where is the time of encounter and where B*R and B*T are projections of 
the targeting vector B along the ordinate and abscissa of the R-T coordinate 
system defined in Appendix F . The vector B is the normal vector -directed 
from the center of the target planet to the approach asymptote. 

Expanding Eq. (2.80) about the nominal trajectory and neglecting terms 
of order higher than the first yields 

T = T * + [w]" • (2 * 82) 

When the deviations Sr and Sv are small enough so that this linearization is 
valid, the deviation of the targeting vector before and after the guidance 
correction can be written 


ST 

= M 1 Sr + 

5v 

(2.83) 

ST + 

= Sr + 

M 2 6v + 

(2.84) 


where 

M 1 = [3T/3?J* 

M 2 = [3T/3v]* 

and ( ) and ( ) indicate quantities just before and after the correction. 
Solving Eqs . (2.83) and (2.84) for the required correction gives 

- + -1 + 

AV = Sv - Sv = M 2 [ST - ST ] (2.85) 

Substituting back in for ST~ gives 

AV = M 2 1 [5T + - (M 1 Sr + M 2 Sv)J . (2.86) 

Assuming a full correction is made, ST + = 0 and Eq. (2.86) becomes 
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AV 



Sr 


I Sv 


Recalling that 


Eq. (2.87) can be written 



AV = Bx 


where B 


r-M ~ 1 m ! -ij. 

I 1 i 


(2.87) 


( 2 . 88 ) 


Updating the State Estimate Error Covariance Matrix, P 

Eq. (2.88) expresses the velocity correction required to null the tar- 
get error as a linear function of the deviation of the actual state from the 
nominal state. However, since the actual state is never known, it is neces- 
sary to compute the commanded guidance correction from the best estimate of 
the state x, i.e. , 


AV = Bx. (2.89) 

Considering execution errors, the actual correction will be 


AV 1 = AV +. n . (2.90) 

where n, the execution error vector, is assumed to be distributed with 
zero mean and known covariance N, i.e., 


N = E[nn T J 


(2.91) 


as defined by Eq. (2.47), 

If a full correction is made, the new state estimate becomes 
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x = x + JAV 

where J = is a 6x3 compatibility matrix, 

mate after the correction is given by 


(2.92) 


The error in the state esti- 


v+ %_ 
x = x + Jn . 


(2.93) 


The state estimate covariance matrix after the correction is defined 


P + - BIST ^ +T J. 


by 


(2.94) 


Substituting for x + yields 


_+ _ r ^- ^-T T T '\»—T T T 

P -E.[x x + x nJ + Jrix + Jnn j j , 


(2.95) 


With the assumption that the execution errors and the state estimate 


errors are uncorrelated (i.e., E{x" n T J = 0), the covariance of x + becomes 


where P = 


E[x x ] and N 


P + = P + jnj t 


E[nn T J. 


(2.96) 


Updating the State Error Covariance Matrix, S 

Consider the update of the state error vector 


x = x + JAV » . 


Using Eqs . (2.89) and (2.90) gives 


x = x + JBx + J n . 


Substituting for x yields 


(2.97) 


(2.98) 


X 


x + JB(x-x) + Jn 


(2.99) 



36 


and thus 


+ 


x 


(I + JB)x - JBx + Jn • 


( 2 . 100 ) 


The state error covariance matrix after the guidance correction can be written 


S 


+ 


= E{x + x +T ] 


( 2 . 101 ) 


Substitution of Eq. (2.100) into Eq. (2.101) yields 

T T 

S + = E [ ( I + JB)x"x" (I + JB) T - JBx" x" (I + JB) T 


T T 

- (I + JB)x x B^ J^ - Jnx (I + JB J 


( 2 . 102 ) 


T T 

T T a— T T T T 

+ JBx x B J - Jnx B J + Jnn J ] 


Before simplifying Eq. (2.102), consider the correlation of the state esti- 

% 

mate error x and the state deviation x, i.e., 


E [x x T ] = E£x (x + x) T ] 

E{x x T ) = E[x x T J t E{x x T ] 


(2.103) 


Using the property that the state estimate and the error in the state esti- 

(14) 

mate are uncorrelated , Eq. (2.103) becomes 

E{x x T ] = P. (2.104) 


Furthermore, the execution error is assumed to be statistically independent 
of both the state deviation and the state estimate error, i.e., 

T 

E{nx ] = E£nx ] = o. (2.105) 

Using Eqs . (2.104) and (2.105) and the symmetry property of the covariance 
T 

matrix (i.e., P = P), Eq. (2.102) can be simplified to yield 
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S + = (I + JB)(S - P )(I + JB) T + P~ + JNJ T . (2.106) 

Recalling Eq. (2.96) 

S + = (I + JB)(S" - P")(I + JB) T + P + . (2.107) 

Eq. (2.107) is the update equation for the state deviation covariance matrix 
when a guidance correction is made. 

Velocity Correction Covariance Matrices 

Although the true velocity correction, AV, is never known its covari- 
ance matrix can be expressed as 

T 

ElAV AV T J = E{Bx"x~ B T J = BS _ B T • (2.108) 

Likewise, the covariance matrix of the error in the estimate of the 
required correction, U, can be expressed by 

U = E [ ( AV - AV) (AV - AV) T J. (2.109) 

Substitution of Eqs . (2.88) and (2.89) into Eq. (2.109) leads to the follow- 
ing development: 

_ rn 

U = E[(Bx - Bx)(Bx - Bx) ] 

U = B E[(x - x)(x - x) T ]B T (2.110) 

U = B P~B T . 

The covariance matrix of the commanded velocity correction, AV, is defined 
by 

V = E [AV AV T J- 


( 2 . 111 ) 
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Substitution of Eq. (2.89) into Eq. (2.111) yields 

V = E[B x x T B T ] 

V = BE[x" - x")(x~ - x~) T ]B T 

V = BE [ x x - x x T - x'x' + x~x _ ]B T 

V = B[S” - P"]B T 

_^_T 

since E[x x ] = 0. 


( 2 . 112 ) 


Statistics of the Magnitude of the Velocity Correction 

Under the assumptions of this study , the commanded velocity correc- 
tion can be characterized by the following statistics : 

E{AV] = <f> (2.113) 

EjAV AV T ] = V (2.114) 

where the expectation is taken over the ensemble of all possible trajectories. 
However, to evaluate guidance propulsion requirements the statistics of the 
magnitude of the correction, |av|, are needed. Unfortunately the mean and va- 

. A 

riance of |AV| are not readily obtained from Eqs. (2.113) and (2.114). Guid- 
ance requirements in Chapter 4 are expressed in terms of the RMS (root -mean- 
square ) value of |AV| given by < ‘ 12 ’ ) 


AV = /Trace V . 
rms 


(2.115) 


Statistics of the Targeting Error 

The actual error in the targeting vector is given by 


<5T = T - T* = Bx. 


(2.116) 
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The covariance of the error is defined by 

cov{6T] = E[6T 6T T J = E[Mx x T M T J (2.117) 

where M = " = [M 1 j 

Since M is deterministic, Eq. (2.117) can be written 

cov [6T J = MSM T . (2.118) 


A convenient parameter for evaluating targeting uncertainty is the target RMS 
defined by 


T 

rras 


= A: 


race 


(MSM ) 


(2.119) 


Guidance Algorithm 

ihe computational algorithm for the guidance scheme derived in the pre- 
vious section can be summarized as follows: 

Given x, P and S at the correction time t, : 

k 

(1) Evaluate X = X* t x, M_ L = QfJ “ and M 2 

(2) Compute T from X and evaluate 6T = T - T*. If | 6T | < e , where e is a 
small tolerance, skip to step (7), otherwise go to step (3). 

(3) Solve for the linear correction 

AV = -M 2 _1 (ST 

(4) Update the state vector X = X . , + JAV 

old 

(5) Evaluate and at the new state X and repeat steps (2) through (5). 

(6) Form B = [-My 1 M j -I] 
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(7) Compute covariance matrix N from Eq. (2.47) 

(8) Update P : P + = P~ + JNJ T 

(9) Update S : S + = [I t JB] [S~ - P~][I + JBJ T + P + 

(10) Compute V and Trace V : 

V = B[s" - P”.]B T 

(11) Compute U and Trace U : 

u = bp" b t 

(12) Compute AV : AV = /Trace V 

rms rms 

(13) Compute T : T = /race MS + M T . 

rms rms 

The derivations presented, in this chapter constitute the theoretical 
basis for the navigation analysis reported in Chapter 4 and the guidance 
analysis presented in Chapter 5. In selecting algorithms and error models 
to simulate the navigation and guidance process, special emphasis was placed 
on methods which are rapid enough to be practical, yet are realistic enough 
to be meaningful. The computational algorithms and models were incorporated 
into two computer programs which are described in the following chapter. 



CHAPTER 3 


COMPUTATIONAL DEVELOPMENT 

The computational algorithms derived in Chapter 2 were incorporated 
into two special purpose computer programs to facilitate the navigation and 
guidance analysis. The programs, designated STEP V and STEP VI, were devel- 
oped to satisfy the requirement for efficient research-oriented mission 
analysis computer programs. They differ only in that STEP V performs in a 
conventional batch mode while STEP VI is used in an interactive mode. Both 
programs were designed for operation on the CDC 6600 computer. This chapter 
contains a general description of the two programs. A more complete documen- 
tation is contained in Reference 15. 

3.1 Description of the Programs 

STEP V is a research-oriented FORTRAN IV computer program designed 
specifically for the preliminary analysis of interplanetary missions. The pro- 
gram simulates the navigation and guidance process for the encounter phase of a 
planetary flyby mission. The nominal and the simulated encounter trajectories 
are generated simultaneously through parallel numerical integrations of the 
equations of motion. A state transition matrix is generated by numerical in- 
tegration of the perturbation equations. All integrations are performed with 
a fourth order Adams-Moulton predictor-corrector scheme using a fourth order 
Runge— Kutta starter. Observational data are simulated at specified intervals 
by adding randomly generated errors to observations computed on the "true” 
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(simulated) trajectory. A linear sequential estimation scheme is used to 
process the observations and to provide an estimate of the state after each 
observation. The error in the estimate is computed by differencing the true 
and the estimated states. 

Guidance corrections , which null the estimated error in the target 
vector, can be executed at specified times. Options are provided to plot 
time histories of errors in the estimate for the simulated case as well as 
pertinent properties of the error covariance matrix associated with the state 
estimate. A basic flow chart of STEP V logic is presented in Figure 7. 

STEP VI, an adaptation of STEP V, can perform all of the computational 
capabilities described above while operating in an interactive mode. The pro- 
gram was designed to take advantage of the interactive graphics capability of 
the Control Data Corporation 252 Display Console at The University of Texas at 
Austin. Primary program functions are shown in Figure 8 and described below. 

Overall program control is from the console keyboard. Executions of 
STEP VI subprograms which control (1) problem input, (2) input data manipula- 
tion, (3) problem setup, (4) problem execution, (5) output data plotting, and 
(6) display storage are accomplished by commands from the console keyboard. 

In addition, certain display instructions can be input directly from the key- 
board to the CDC 252 system. Since the program normally occupies central 
memory for the entire working session, storage requirements were minimized by 
dividing the program into three functional segments. The control program, 
which converts the keyboard commands into STEP VI subroutine calls , is the 
main segment and is always resident in central memory. Segment A, comprised 
of subroutines which perform functions 1, 2 and 5 shown in Figure 8, is 
loaded from disk when one of its subroutines is first called by the main seg- 
ment. Segment B replaces Segment A in central memory when a subroutine which 
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accomplishes functions 3 or 4 is called by the main segment. To conserve 
storage, output data generated by Segment B is stored on disk and read off 
by Segment A when the plotting command is given. 

3.2 Description of Interactive Program Operation 

An interactive session is initiated by reading the program into the 
CDC 6600 system input queue as a batch job. When the program is assigned a 
system control point and update compilations performed (if any), the main 
segment is loaded and a message is displayed on the CRT to indicate that the 
program is waiting for a keyboard command. The mission analyst causes the 
first problem to be loaded by keying in a simple command (the letter "I”). 

The command causes Segment A to be loaded, the input routine to read the first 
problem, and a problem identification message to be displayed. The analyst 
may procede with the first problem or he may sequence to any other problem by 
repeating the input command. When the desired problem has been read, the 
mission analyst can display all input parameters and execute changes to them 
from the keyboard. Changes are accomplished by typing in two identifying num- 
bers and the new value of the parameter. When a change is executed, the new 
value is displayed to the right of the original input value. When the analyst 
is satisfied with all input parameters, he keys in the setup command which 
causes Segment B to be loaded, all problem setup operations to be performed, and 
the message "GO" to be displayed. 

Numerical integration of the trajectories (nominal and simulated), and 
thus execution of the problem, are controlled by the keyboard command "GO, II, 
12." The parameter II is an integer which specifies the number of integration 
steps to be performed. The parameter 12 is an integration status display indi- 
cator. If the status display option is selected, the current time-to-go (to 
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encounter) and distance-to-go are displayed on the CRT during execution. If 
a guidance correction time has been specified, problem execution will stop 
when the correction time is reached. Pertinent guidance data are automati- 
cally displayed at that point. The analyst then has the option to execute 
the correction and continue, or skip the correction and return problem control 
to the keyboard. When the specified integration is complete, problem control 
is automatically returned to the keyboard. 

At this point the analyst may plot the data generated on the run just 
completed by giving a simple plot command. The plot command calls Segment A 
back into memory and causes selected data to be automatically scaled and plot- 
ted on the CRT. By keying in the appropriate commands the plotted data can be 
recorded directly on microfilm and/or stored on magnetic disk for recall later 
in the session. Each plot is tagged with the time, the date, and a unique 
identification number. The tag information is also recorded on the printed 
output . 

After completion of a problem the mission analyst may terminate the 
session, rerun the problem with new input, or read in a new problem. If a 
Monte Carlo analysis has been specified in the input data, the analyst need 
only key in the setup command to generate a new set of initial state-errors. 

Central processor time for a typical planetary approach problem aver- 
ages less than one minute. However, since the program operates in a multi- 
programming environment, the elapsed session time for the same problem will 
vary from 2 to 5 minutes. Maximum storage required is 70,000 octal locations. 



CHAPTER 4 


APPROACH NAVIGATION ANALYSIS 

A comprehenisve analysis of the approach phase of Jupiter, Saturn 
and Uranus encounter was made to define reference approach navigation accu- 
racies at those planets. For the purposes of the analysis, the approach 
phase is defined as that segment of the encounter trajectory extending from 
the SOI (sphere of influence) to periapsis. This is the critical period in 
the mission profile when the state of the spacecraft relative to the planet 
must be determined with sufficient accuracy to permit the execution of one 
or more guidance corrections . 

When the spacecraft reaches the target planets' SOI, the uncer- 
tainty in its heliocentric state will have been improved by many months of 
Earth-based radar tracking data. However, the state of the spacecraft rela- 
tive to the target planet will be known to a much lesser degree of accuracy 
due to the relatively large uncertainties remaining in the ephemeris of the 
outer planets . Approach navigation must reduce this uncertainty to accep- 
table levels before an effective guidance correction can be performed. Nav- 
igation errors which remain at the time of a guidance correction will result 
in an incorrect guidance execution. If the guidance errors remain uncor- 
rected, they will propogate through the encounter resulting in much larger 
post-encounter errors. 

4.1 Assumed Error Sources 
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During the planetary approach phase , navigation errors result from 
errors in the observations as well as from errors in the assumed dynamical 
model. For this study, the primary source of error was assumed to be a 
random noise in the observation. Other sources of observation error not con- 
sidered in this analysis include errors in the Earth-based tracking station 
location and center-finding errors inherent in the planet and sun sensors. 
Since these errors are systematic in nature, they can be reduced signifi- 
cantly by including the biases as parameters to be estimated. The principal 
source of error in the dynamical model was assumed to be the uncertainty in 
the position of the outer planets . Another source of model error not con- 
sidered here is due to the uncertainty in the mass of the planet. As with 
the other systematic errors the uncertainty can be reduced by including the 
planet mass as a parameter to be estimated. 

Nominal one-sigma values for the assumed error sources are listed 

in Table 1. These data are based on projected improvements in the accuracy 

( 16 ) 

of the Deep Space Network . Onboard angle errors are considered to be 
typical of optical measurement systems currently under development 

at the Jet Propulsion Laboratory . Planet position errors are based on an 
expected 0.2 arc-sec uncertainty for each planet^ ^ . All error sources used 
in this study are considered to be on the conservative side of projected 
values for the late 1970s. 

4.2 Reference Navigation Accuracy 

Computer simulations of the navigation process for Jupiter, Saturn 
and Uranus encounter were made for the nominal 1977 Grand Tour mission de- 
rined in Chapter 1 and Appendix A. For the purpose of defining reference 
navigation accuracy data, the following ground rules were adopted: 



Table 1 Nominal One-Sigma Values for Assumed Error Sources 


Source One-Sigma Values 

Observations : 

Range 

Range-rate 
Sun-planet angle 
Star-planet angle 

Planet Position: 

Jupiter 750 km 

Saturn 1400 km 


15 m 

.5 mm/sec 
10 arc-sec 
10 arc-sec 


Uranus 


2700 km 
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(1) Continuous Earth-based tracking data are available from a 
single station except when the spacecraft is occulted by the 
Earth . 

(2) An equivalent data rate of one observation every 144 minutes 
(0.1 days) for all observation types. A realistic Earth-based 
data rate (i.e., one observation per minute) is simulated by 
scaling the nominal observation variance by the ratio of the 
realistic data rate to the assumed data rate (i.e., the ratio 
1/144). 

(3) The error in the location of the Earth-based tracking station 
is zero . 

(4) The initial state error is zero (initial values of the true 
and nominal state deviation vectors are equal). 

(5) The apriori estimate of the state deviation is zero (x =0). 

o 

The apriori covariance matrix of the estimate, P , is a diag- 
onal matrix defined in Table 2. 

Reference navigation accuracy data are presented in Figures 9 through 
27. These results are based on, and subject to the limitations of, the linear 
estimation theory derived in Chapter 2. The data presented are of three basic 
types: (1) state estimate uncertainty derived from the state estimate error 

covariance matrix, (2) target error ellipses derived from the target estimate 
covariance matrix, and (3) target error data obtained from mission simula- 
tions . 

State Estimate Uncertaint y 

Reference state estimate uncertainty during Jupiter, Saturn and 
Uranus approach are presented in Figures 9 through 11. These data consist 
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Table 2 Diagonal Elements of the Apriori State Estimate 
Covariance Matrix P 

o 


Element 

Jupiter 

Saturn 

Uranus 

P (km 2 ) 

°11 

5.63 (5) 

2.00 (6) 

7.62 (6) 

P (km 2 ) 

°22 

5.63 (3) 

2.00 (6) 

7.62 (6) 

P (km 2 ) 

°33 

5.63 (5) 

2.00 (6) 

7.62 (6) 

P (km 2 /sec 2 ) 

44 

4.00 (-8) 

4.00 (-8) 

8.00 (-8) 

P (km 2 /sec 2 ) 

°55 

4.00 (-8) 

4.00 (-8) 

8.00 (-8) 

P (km 2 /sec 2 ) 

66 

4.00 (-8) 

4.00 (-8) 

8.00 (-8) 

P o 7? (km 2 ) 

5.63 (5) 

1.96 (6) 

7.29 (6) 

P (km 2 ) 

°88 

5.63 (5) 

1.96 (6) 

7.29 (6) 

P (km 2 ) 

® nn 

5.63 (5) 

1.96 (6) 

7.29 (6) 



of RMS (root-mena-square) values of the uncertainty in the space-craft's 
planet ocentric position and velocity and the heliocentric position of the 
planet, as well as an overall measure of state uncertainty. RMS uncer- 
tainty in position, velocity and planet position are defined by the square 
root of the sum of the component variances, i.e.. 


Position RMS 

= 

( p n 

+ 

P 22 

+ 

p )”* 
33' 

Velocity RMS 

= 

( p 44 

+ 

P 55 

+ 

P V 5 
66' 







Jg 

Planet Position RMS 

= 

(P?? 

+ 

CO 

CO 

p-i 

+ 

P ) 
99' 


(4.1) 


where the variances P^, ^22 , *’ , ’^99 are *^ e ^ a § ona -L e l eTnen ts ° f the s t a te 

estimate error covariance matrix. Overall state uncertainty is expressed as 

the "entropy" of the system. The concept of entropy as a measure of state 

(17 ) 

uncertainty is well accepted in the field of information theory . In 

f 1 8 ) 

1967 Potter and Fraser ' derived the following expression for the entropy 


of a sequential estimator as a function of the determinant of the state es- 


timate covariance matrix: 


Entropy = |p| + E q (4.2) 

where E = (n/2)loe„ (2ire), n is the dimension of the state vector, and e is 
o 2 

the base of the natural logarithm. Although entropy as used here has no 
physical significance, it is considered a convenient figure-of-merit for 
overall performance of the navigation process . 

State Estimate Uncertainty During Jupiter Approach 

Reference navigation accuracy data for the state estimate during 
Jupiter approach are presented in Figure 9 for (1) range rate, RR, (2) on- 
board angles, OBA, and (3) range-rate plus onboard angles, RR + OBA. 
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Time histories of position RMS during Jupiter approach are shown in 
Figure 9a. From Figure 9a it can be seen that the OBA information alone is 
not sufficient to reduce significantly the uncertainty in spacecraft position. 
In fact , the information rate from onboard angles is not adequate to prevent 
a growth in the position RMS during the period from E-30 (encounter minus 30 
days) to E-6 days. Range-rate data alone give much better position accuracy 
for the same period. A significant overall improvement in position accuracy 
is achieved by combining range-rate' observations with the onboard angle ob- 
servations . 

Velocity RMS time histories for Jupiter approach are shown in Fig. 

9b. These data are characterized by a sharp peak in the velocity uncertainty 
just prior to encounter, followed by a rapid decrease at encounter. This 
phenomena was observed to occur at each planet for all observation types, and 
can be explained by the following considerations. The change in the covari- 
ance matrix P is governed by two effects, viz., a decrease caused by the in- 
formation content of the observations and an increase due to the dynamical 
propagation of the covariance matrix. For a short period just prior to en- 
counter the information rate is not great enough to offset the growth in the 
velocity elements of the covariance matrix, thus causing the velocity RMS to 
increase sharply before the information term begins to dominate again. The 
onboard angle data reduce the uncertainty slightly up to E-12 days at which 
point the characteristic velocity uncertainty growth begins. Range-rate 
data give good velocity information up to E-8 days before the velocity un- 
certainty growth takes over. Onboard angle data combined with range-rate 
data give only slight improvement along most of the curve but virtually 
eliminate the velocity uncertainty peak . 

Time histories of Jupiter position uncertainty are presented in Fig. 
9c. These data show that range-rate observations alone provide a significant 
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(9a) 


(9b) 


(9c ) 


TIME FROM ENCOUNTER (Days) 




Figure 9 State Estimate Uncertainty During Jupiter Approach 
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improvement in the position uncertainty of Jupiter while onboard angle mea- 
surements alone provide no information. However, range-rate and onboard 
angle measurements compliment each other to provide a significant overall 
improvement in Jupiter position uncertainty over the range-rate only case. 

State Estimate Uncertainty During Saturn Approach 

Reference navigation accuracy data for the state estimate during 
Saturn approach are presented in Fig. 10 for (1) range-rate, RR (2) onboard 
angles, OBA , and (3) range-rate plus onboard angles, RR + OBA. 

Position RMS time histories. Fig. 10a, show that the onboard angles 
data give slightly better position accuracy than do the range-rate data up 
to about E-15. After E-15 the range-rate data is superior and yields much 
better results near periapsis. As was the case at Jupiter, a combination of 
range-rate and onboard angle measurements gives much better position accu- 
racy than either range-rate or onboard angles alone. Overall position ac- 
curacy is considerably better at Jupiter approach than at Saturn approach. 

Velocity RMS time histories for Saturn approach are presented in 
Figure 10b. As can be seen from Figure 10b range-rate data give a much bet- 
ter velocity estimate than the onboard angle data. For onboard angle data 
the velocity uncertainty growth starts at about E-5, reaches a peak of 125 
km/day and then drops sharply at periapsis. The combined range-rate and on- 
board angle data give a steady and rapid decrease in velocity uncertainty. 

As can be seen from the time histories of Saturn position uncer- 
tainty in Fig. 10c, both range-rate data and range-rate plus onboard angle 
data provide an improved estimate of planet position during approach. The 
onboard angles alone provide virtually no information regarding planet po- 


sition . 
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State Estimate Uncertainty During Uranus Approach 

Reference navigation accuracy data for the state estimate during 
Uranus approach are presented in Fig. 11 for (1) range-rate, RR, (2) on- 
board angles, OBA, and (3) range-rate plus onboard angels, RR + OBA. 

From the position RMS time histories, Fig. 11a, it can be seen that 
the onboard angle data provide better position accuracy at Saturn than do 
the range-rate data up toE-10. After E-10 , the range-rate data becomes supe- 
rior. Range-rate plus onboard data provide a significant improvement over 
both the range-rate case and the onboard angle case after E-25. 

Velocity RMS time histories for Uranus approach are presented in 
Fig. lib. As was the case with Jupiter and Saturn, range-rate data provide 
better velocity accuracy during Uranus approach than do the onboard angle 
data. Unlike Jupiter and Saturn, however, combining onboard angle data with 
range-rate data results in a significant decrease in velocity uncertainty 
throughout Uranus approach. 

Time histories of Uranus position RMS are presented in Fig. lie. 

Only slight improvement in Uranus position uncertainty can be achieved prior 
to periapsis passage with range-rate data. Although the onboard angle data 
alone provide virtually no Uranus position information, they significantly 
reduce the uncertainty when combined with range-rate data. 

Target Estimate Accuracy 

Reference navigation accuracy data for Jupiter, Saturn, and Uranus 
approach target parameters B-R and B-T are shown in Figures 12 through 17. 

The data consist of target error ellipses in B-R - B-T target space. The 
error covariance matrix associated with the estimate of the target parame- 
ters is used to generate a contour of constant probability which, for the 
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linearized problem considered here, is an ellipse centered at the estimated 
values of B.R and B-T (see Appendix E). The reference ellipses presented in 
this chapter are based on one-sigma values for B.R and B-T. Correspondingly, 
there is a 40 percent probability that the actual (true) values of the tar- 
geting parameters lie within the error ellipse. Error ellipses corresponding 
to any probability value can be derived from the reference ellipse by 
multiplying the ordinate and abscissa by the factor [2 In(P^-l)]. For 86 and 
99 percent ellipses this quantity has the values of 2 and 3 respectively. 

Jupiter Target Error Ellipses 

Reference navigation accuracy data for Jupiter approach target pa- 
rameters B-R and B-T are shown in Figures 12 and 13. In Fig. 12 the stan- 
dard target error ellipses are compared at E-30, E-15, and E-5 days. These 
times were selected as representative of the period during which an approach 
guidance maneuver will be considered. Fig. 13 demonstrates the relative 
navigation accuracy of various combinations of the observation types. The 
error ellipses in these figures represent the dispersion of the true target 
points (B-R and B-T) relative to the estimated target point. (Recall that 
the target point is the point where the approach asymptote pierces the tar- 
get plane. See Figure 39 for a definition of error ellipse nomenclature.) 

Target error ellipses for range-rate, RR, observations only are 
shown in Figure 12a. These data indicate a slight decrease (100-200km) in 
both coordinates from E-30 to E-15 days. During the next ten days of track- 
ing the B*T coordinate decreases an order of magnitude while the B-R coor- 
dinate remains virtually unchanged. Since both the spacecraft approach 
asymptote and the orbit plane (and thus T) lie very close to the ecliptic 
plane at Jupiter, most of the range-rate information goes into the estimate 














of the B-T component. B-R is essentially an out-of-plane component and is 
therefore more difficult to estimate with range-rate data alone. 

In Fig. 12b the orbit determination ellipse for the onboard angles 
alone (sun-planet angle, a, and star-planet angle, 3) are presented at E-30, 
E-15, and E-5 days. Note that onboard angle data give about a 370 km reduc- 
tion in both B*R and B-T coordinates from E-30 to E-15. From E-15 to E-5 
the B-T (in-plane) information is slightly better than the B-R (out-of-plane) 
information. 

In Fig. 12c error ellipses for range-rate and onboard angle data 
show a rapid reduction in both B-R and B-T components when the data types 
are combined. Navigation accuracy is still better in the B-T direction than 
in the B-R direction due to high quality range-rate information in the in- 
plane direction. 

Target error ellipses for various combinations of observation types 
are compared in Fig. 13. In addition to the data previously presented in 
Fig. 12 the data in Fig. 13 include error ellipses for range-rate plus sun- 
planet angle and range-rate plus star-planet angle. The sun-planet angle, 
a, being basically an in-plane angle, does not reduce the B-R component sig- 
nificantly over the range-rate only case. However, it does provide better 
in-plane information than range-rate alone as evidenced by the reduction in 
B-T at E-30 and E-15. At E-5 the difference between the range-rate plus sun- 
planet angle case and the range-rate case is negligible. 

Saturn Target Error Ellipses 

Reference navigation accuracy data for Saturn approach target pa- 
rameters B-R and B-T are shown in Figures 14 and 15. In Fig. 14 the stan- 
dard target error ellipses are compared at E-20, E-10 and E-2 days. The 
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relative accuracy of various observation combinations are presented in Fig. 
15. 

Target error ellipses for range-rate observations are shown in Fig. 
14a. From that figure, it can be seen that virtually no target information 
is obtained from range-rate observations from E-20 to E-10 days. During the 
next 8 days of tracking a significant reduction in B*T is achieved with no 
significant decrease in B*R uncertainty. In contrast to the range-rate data, 
onboard angle data provide good target information throughout Saturn approach 
as can be seen from Fig. 14b. The error ellipses at E-20 and E-10 are nearly 
circular indicating that the onboard angles provide information of almost 
equal value in both the in-plane and out-of-plane directions . Range-rate 
plus onboard angle data provide about the same B*R accuracy as the onboard 
angle only case, but much improved B*T accuracy. 

In Figure 15 the relative accuracy of the various data combinations 
is compared at E-20, E-10, and E-2 days. From Fig. 15 it can be seen that 
the star-planet angle, 8, when combined with range-rate, p, does not reduce 
the B*R component at Saturn significantly up through E-10 days. At E-2 
days, however, the combination p + 3 gives a target accuracy comparable to 
the p + a + 8 case and better than the a + 8 case. As was true at Jupiter 
approach, range-rate data alone provides a poor estimate for B»R. 

Uranus Target Error Ellipses 

Reference navigation accuracy data for Uranus approach target pa- 
rameters B*R and B*T are presented in Figures 16 and 17. In Fig. 16 the 
standard target error ellipses are compared at E-20, E-10, and E-l days. 

The relative accuracy of various observation combinations are presented in 


Fig. 17. 




( a ) RR (b) OBA (c) RR+OBA 


Figure 16 Navigation Target Error Ellipses at Uranus Approach (RR, OBA, RR + OBA) 
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Target error data for range-rate observations are shown in Fig. 16a. 
Range-rate data provides practically no target information from E-20 to E-10 
days. The error ellipses at E-20 and E-10 are very nearly circular with a 
radius of about 2000 km. From E-10 to E-l the B-T uncertainty is reduced to 
below 150 km with essentially no decrease in B-R uncertainty. Onboard angle 
data, however, give a uniform decrease in both B-R and B-T as shown in Fig. 
16b. Combining range-rate data with onboard angle data results in a smaller 
B-T uncertainty than with onboard angle data only (Fig. 16c). 

In Fig. 17 the relative accuracies of various data combinations are 
compared at E-20, E-10, and E-l days. Significantly, the range-rate data 
are ineffective in reducing target errors until well into the approach phase. 
Adding the star-planet angle to the range-rate measurement (i.e., p + 6) de- 
creases B-R uncertainty significantly while leaving B-T uncertainty virtually 
unchanged. Similarly, adding the sun-planet angle to the range-rate measure- 
ment (i.e., p + a) results in a large decrease in B-T uncertainty while 
leaving B-R uncertainty virtually unchanged. A combination of all three ob- 
servations, p + a + B, yields good target accuracy in both the B-R and B-T 
components . 

Tables 3, 4 and 5 summarize in tabular form the navigation data pre- 
sented in the previous sections. 

Overall State Uncertainty 

Overall state uncertainties , as measured by the entropy of the esti- 
mation process, are presented in Fig. 18 for Jupiter, Saturn and Uranus ap- 
proach. These data reflect the fact that although onboard angle measurements 
alone contain limited overall state information, they do provide valuable 
complimentary information when combined with range-rate data, thus making pos- 
sible significantly improved navigation accuracy. 



TARGET ERROR VARIANCES: 

°B.R (km2) 
al m (km 2 ) 


a £ ( days 2 ) 


5.2 (+5) 4.2 ( + 5 ) 2.3 ( + 5) 3.K+5) 7.9(+4) 1.6 (+4) 3.1(+5) 8.3(+3) 3.0(+3) 

3 . 4 (+5 ) 4.2 (+5 ) 3.6(+4) 2.l(+5) 7 . 8 ( +4 ) 7.8(+3) 4.4(+2) 1.7(+4) 3.7(+2) 

1 . 2 (-7) 6 . 7 (-7) 2 . 8 ( -8 ) 7 . 3 (-7 ) 6.o(-6) 1 . 6 (- 8 ) 6.9(-9) l.K-6) 9.2(-ll) 












a le 4 Comparison of Navigation Accuracies at Saturn Approach 


Observation “ 

STATE UNCERTAINTY: 


E-20 

"OBA |~RR + OBA 


E-10 

"OBA 


RR + OBA" 


Position RMS (km) 

1479 

1450 

315 

1459 

1466 

147 

1009 

1440 

51 

Velocity RMS (km/day) 

17 

29 

16 

13 

25 

11 

13 

36 

6 

Planet Position RMS (km) 

1957 

2425 

1470 

1810 

2423 

1206 

1209 

2422 

389 

Entropy 

60 

88 

55 

51 

85 

44 

37 

80 

35 
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labia 5 Comparison of Navigation Accuracies at Uranus Approach 


OBSERVATIONS 


OBA RR + OBA RR 


OBA RR + OBA RR 


OBA RR + OB' 


STATE UNCERTAINTIES: 


Position RMS (km) 
Velocity RMS (km/day) 
Planet Position RMS (km) 
Entropy 


TARGET ERROR ELLIPSE: 


a (km) 
b (km) 
(deg) 


IARGET ERROR VARIANCES: 

° B*R (km2) 

oi „ (km 2 ) 


o 2 ( days 2 ) 


4.1 (+6) 3. 7 ( + 5) 3.5 ( + 5) 3.9(+6) 8.4( + 4) 6.8(+4) 3.7( + 6) 1.3( + 3) 1.3( + 3) 

4.0 (+6) 3.7 (+5 ) 5.8 (+4) 4.0(+6) 8.4(+4) 1.3 (+4) 1.0(+5) 3.3(+3) 7.7(+2) 

2.8 (-4) 1.1 (-4) 4.2 (-5) 9.3(-6) 2.3(-6) 1.4 (-7) l.l(-9) 1.6(-6) 2. 5( -12 




Encounter 







73 


Target Errors for Simulated Encounter 

The two types of data presented previously are based on assumed sta- 
tistical models for the error sources and represent predicted navigation er- 
ror statistics for state and targeting components. The actual navigation 
errors for any single mission will depend upon observation errors which ac- 
tually occur during the course of that single mission. Although the ac- 
tual navigation errors are never known apriori, simulation of the orbit de- 
termination process using assumed error sources provides a means of studying 
the estimation process and evaluating the validity of the predicted error 
statistics. Example time histories of actual target errors for simulated 
encounters with Jupiter, Saturn and Uranus are shown in Figs. 19 through 27. 
For these simulations both the initial state error and the initial state 
estimate are, zero, i.e., x q = x q = 0. The target error is defined as the 
true (simulated) value minus the estimated value of the components of the 
target vector. 

Time histories of the target vector at Jupiter for the simulated 
mission are presented in Figs. 19, 20 and 21 for RR, OBA and RR + OBA cases, 
respectively. For all three cases the estimation process converged in a 
manner consistent with the predicted navigation accuracies listed in Table 3. 
For the RR case the characteristic velocity uncertainty growth. Fig. 9b, is 
reflected in an increased B*R error just prior to encounter, Fig. 19a. The 
velocity uncertainty for the OBA case. Fig. 10b, is reflected in increased t 

e 

error near encounter. Fig. 20c. A combination of RR and OBA, Fig. 19c, pro- 
vides a significant improvement in estimation accuracy over the RR and OBA 
cases. 

Figs, 22, 23 and 24 consist of time histories of the target errors 
at Saturn for the RR, OBA and RR + OBA cases respectively. For all three 
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cases the error histories are consistent with the predicted navigation accu- 
racies listed in Table 4. The RR data. Fig. 22, exhibit a B*R error growth 
just prior to encounter which coincides with the increase in velocity error 
covariance shown in Fig. 10. For the OBA case the velocity uncertainty is 
reflected in the error in time of encounter, t , Fig. 23c. Combining RR and 
OBA, Fig. 24, results in a significant improvement in all three components 
of the target vector. 

Time histories of target errors at Uranus are presented in Figs. 

25, 26 and 27 for the RR, OBA and RR + OBA respectively. From Fig. 25a it 
can be seen that the errors in B-R and B*T are not reduced significantly un- 
til encounter, a result which is consistent with predicted navigation accu- 
racies listed in Table 5. The time-of-encounter estimate is very good after 
about E-10 days. The OBA data. Fig. 26, provide a better estimate for B*R 
and B'T but give a poor estimate for t just prior to encounter. This re- 
sult is consistent with the predicted navigation accuracies presented in 
Table 5. A combination of RR and OBA data. Fig. 27, give a significant im- 
provement over both RR and OBA data. 




Encounter I \ Encounter I I Encounter 
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TIME TO ENCOUNTER (Days) 



Figure 20 Time Histories of Target Errors at Jupiter Approach (OBA) 
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( 21a ) 


( 21b ) 


( 21c ) 


TIME TO ENCOUNTER ( Days ) 



Figure 21 Time Histories of Target Errors at Jupiter Approach (RR + OBA) 
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(23a ) 


(23b) 


(23c ) 


Figure 23 


TIME TO ENCOUNTER ( Days ) 



TIME FROM SOI ( Days ) 

Time Histories of Target Errors at Saturn Approach (OBA) 


















CHAPTER 5 


APPROACH GUIDANCE ANALYSIS 


Because of the extreme sensitivity of the Grand Tour mission to er- 
rors in the encounter trajectory at Jupiter, Saturn and Uranus, the success 
of the mission will depend to a large extent on the effectiveness of the ap- 
proach guidance procedure. This chapter contains the results of an analysis 
of the approach guidance process for the 1977 Grand Tour Interior mission. 

The approach phase of the encounter at Jupiter, Saturn and Uranus was simu- 
lated to survey and define guidance requirements and target accuracies at 
those planets. A sensitivity analysis was conducted to determine the effects 
of observation accuracy, guidance execution accuracy, and initial state er- 
rors on the approach guidance requirements. Finally a Monte Carlo simulation 
of the approach navigation and guidance process at Saturn encounter was con- 
ducted for the purpose of testing the adequacy of the linear theory, 

5.1 Assumed Execution Errors 

In addition to the assumed error sources listed in Table 1, the 
analysis of this chapter considers errors in the execution of the guidance 
correction. (The execution error model is described in Chapter 2.) The un- 
certainty in the velocity correction is assumed to originate from two inde- 
pendent sources: (l) uncertainty in the magnitude of the correction due to 

a random error in the total impulse applied by the rocket engine, and (2) 
uncertainty in the orientation of the correction vector due to random errors 


84 



in the direction of the thrust vector . Nominal values of one-sigma guidance 
execution errors were assumed as follows : 

Correction magnitude: k = 0.001 (0.1 percent) 

Orientation Angle: y = 0.001 degrees (3.6 arc sec). 

The factor k relates the execution error to the correction magnitude (error = 

K l AV ! ) ’ an<a " tne an gl e Y .is the angle between the commanded correction vector 
and the actual correction . vector . 

5.2 Preliminary Considerations 

A graphical representation of various covariance matrices associated 
with the navigation and guidance process is presented in Figure 28. The ini- 
tial state error covariance matrix, S q , and the initial state estimate error 
covariance matrix, P q , are assumed equal at the sphere of influence (SOI). 

As the orbit determination process procedes , the estimate of the state im- 
proves and the P matrix decreases. The rate of decrease is dependent upon 
the frequency and the information content of the observations. The S matrix 
increases with time as the initial state uncertainty is propagated along the 
nominax trajectory. Since the initial state errors are assumed to be dis- 
tributed about a zero mean, the probability ellipse representing S is centered 
on tne nominal state, X*. The probability ellipse representing P is centered 
at the current state estimate, X. The probability ellipse describing the dis- 
persion of the state deviation estimate, x, defined by the covariance matrix 
S-P, is centered on X-'. When a guidance correction is executed, the velocity 
components of the P matrix are increased due to the uncertainty in the guid- 
ance execution (see Equation (2.96)). The S matrix decreases significantly 
at execution due to the input of navigation information. This decrease is 



X <x lx 


X = ACTUAL STATE 
X = ESTIMATE OF Tl 
X* = NOMINAL STA‘ 
= X - X* 

= X - X* 

_ A 

- X - X 

P = E [xx T ] 

S = E [x x T ] 

S-P = E [x x T ] 
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offset slightly by the added uncertainty resulting from guidance execution 
errors (see Equation (2.107)). 

As derived in Chapter 2, AV and 6T are given bv 

rms rms ° 7 


{Trace [B(S - P^ 1 ]} 3 * (5.1) 

(Trace [M S + (5.2) 

where 

M = £M ± j M 2 ] = [3 T/3 r | 3T/3v] 

B = [-M 2 _1 M j _ | -IJ 
S = E[xxh 
P = E[xx T J 


AV 

rms 


T 

rms 


and where the subscripts ( ) + and ( ) indicate values of the quantities be- 
fore and a±ter the correction. Since the first two components of the target 
vector T (i.e., B-R and B*T) are of primary interest to this analysis, it is 
convenient to adopt the notation T* for the abbreviated target vector, i.e.. 


T T 


B.R 

B-T • 


The state estimate error covariance matrix P and the state error co- 
variance matrix S were assumed to be equal at the initial time (i.e. , at 
entry into the target planet T s SOI). Nominal values of the diagonal elements 
(variances) of the initial state error covariance matrix, P , are listed in 
Table 2. Off-diagonal elements of P^ were assumed to be zero. 

The nominal heliocentric velocity V * of the spacecraft at the out- 
bound SOI is given by 


V 


V* + V * 
P 


(5.3) 
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the nominal heliocentric velocity of the planet. Neglecting the small error 
in the planet T s heliocentric velocity and assuming that the departure correc- 
tion required to account for small heliocentric position errors are negligible, 
the velocity correction required at departure is 


AV 3 = V - V* (5.4) 

where V is the actual planetocentric velocity of the spacecraft. At the SOI 
Eq. (5.4) can be closely approximated by 


AV 


3 



(5.5) 


where V is the hyperbolic excess velocity. Furthermore, at the outbound 
°°D 

SOI, the covariance matrix P becomes small relative to the covariance matrix 
S, thus 


S(t 3 ) - P(t 3 ) s S(t 3 ). 


(5.6) 


The assumptions and approximations just described allow for the computation 
of the covariance matrix associated with the departure correction AV^ with- 
out simulating the mission past the time of the second correction, i.e.. 


cov [AV ] = C {cov [6T ! ] }C 
d t 2 


(5.7) 


where C = 


3 V 


3 T 


is the matrix of partial derivatives required to map from 




approach target space to departure target space (evaluated on the nominal 
path at time t^). Analytical expressions for the matrix C are derived in 
Appendix G. The RMS value of the departure correction is defined by 


AV 

rmS 3 


{Trace 


(cov £ AV g J )}t 


(5.8) 
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5.3 Approach Guidance Requirements and Target Accuracies 

Guidance execution times were surveyed to determine the variation of 
guidance requirements ( AV rms ) and target accuracy ( 6T ^ ms ) during approach to 
Jupiter, Saturn and Uranus. Both a one and a two guidance correction strat- 
e gv was considered for the range -rate (RR) and the range-rate plus onboard 
angle (RR + OBA) cases. 

One Correction Survey 

The variations of AV and 6T T with execution time are shown in 

rms rms 

Figures 29, 30 and 31 for one correction at approach to Jupiter, Saturn 
and Uranus respectively. Since the covariance matrices associated with er- 
rors in the state and the state estimate are assumed to be identically equal 
at the SOI (i.e., P = S ), the AV is zero initially (see Equat ion ( 5 . 1 ) ) . 
The curves increase slowly at first and then increase asymptotically 

as periapsis is approached. The AV^^ curves for the two observation cases 
are very nearly the same with the RR curve falling slightly under the RR + 

OBA curve. The post-correction target accuracy, 6T^ , decreases from SOI 

to encounter with the RR curve substantially above the RR + OBA curve for 
each planet. 

Two Correction Survey 

A survey of two approach guidance corrections was made to determine 
the effect of execution times on two selected indices of guidance perfor- 
mance. The navigation and guidance process was simulated for a matrix of 
five values of first correction time (t ) and six values of second correc- 
tion time These data include the one-correction data previously pre- 

sented since at the SOI the required velocity correction is zero. Repre- 
sentative contours of the two performance indices listed below were drawn in 
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TIME TO ENCOUNTER (Daysl 
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t l ~ t 2 s P ace: 


I 


1 


AV + AV 
rms ^ rms ^ 


T = AV + AV + AV 
2 rms rms ^ rms 3 


AV and AV are the RMS values of the guidance corrections at t, and 
rms. rms 0 1 

1 2 

t respectively. A ^ rms the correction required to null the velocity er- 

3 

ror at departure resulting from target errors remaining after the second cor- 
rection. The departure correction requirement was defined by the correction 
required at outbound SOI to achieve the nominal hyperbolic excess velocity 

V . 

oo 

D 

Surveys of guidance performance index I for the RR case and RR + 

OBA case are presented in Figures 32, 33 and 34 for Jupiter, Saturn and 
Uranus, respectively. These data exhibit a smooth gradient of increasing 1^ 
values as both the first and second guidance times increase. The index is 
much more sensitive to the second correction time than to the first correc- 
tion time. The RR cases have slightly higher I values than do the corres- 
ponding RR + OBA cases. A comparison of the Jupiter, Saturn and Uranus data. 
Figures 32, 33 and 34 respectively, shows that both I and the gradient of 1^ 
increase with each successive planet. For all three planets the correction 
schedule which yields the lowest values of ^ is a single correction at the 
latest time. However, since I does not include the effect of target error 
remaining after the last approach correction, it is of limited interest as a 
performance index. 

A more meaningful survey parameter is the performance index 1 ^. Re- 
sults of the I survey are presented in Figures 35, 36 and 37 for Jupiter, 
Saturn and Uranus respectively. For the RR + OBA case at Jupiter, Figure 
35a, the contour pattern is very similar to the one for 1^, the index I being 



E-15 E-20 E-25 E-30 E-35 E-39.5 

FIRST CORRECTION TIME, tj (Days) 


Figure 32 Survey of Guidance Performance Index Ij, Jupiter 




SECOND CORRECTION TIME, t 2 (Days) SECOND CORRECTION TIME, t 2 (Days) 
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( 33a ) RANGE-RATE PLUS ONBOARD ANGLES 




Figure 33 Survey of Guidance Performance Index 1^, Saturn 
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(34a ) RANGE-RATE PLUS ONBOARD ANGLES 



( 34b ) RANGE - RATE 



Figure 34 Survey of Guidance Performance Index Ij, Uranus 
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Figure 35 Survey of Guidance Performance Index I 2 , Jupiter 
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( 36a ) RANGE-RATE PLUS ONBOARD ANGLES 
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( 36b ) RANGE-RATE 



Figure 36 Survey of Guidance Performance Index I 2 , Saturn 
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( 37a ) RANGE-RATE PLUS ONBOARD ANGLES 



( 37b ) RANGE-RATE 



Figure 37 Survey of Guidance Performance Index I 2 , Uranus 
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much more sensitive to the second correction time than to the first. For a 
given first correction the second correction time for minimum appears to 
be near E-10 days. For the RR case at Jupiter, Figure 35b, a valley centered 
at about E-9 stretches across the survey area. The values for the RR case 

are larger than for the RR t OBA case and exhibit a much stronger gradient. 

The values for the RR case are on the order of two times greater than the 
corresponding I values of the RR + OBA case. 

The survey data for Saturn approach clearly define the region of min- 
imum I^. For the RR + OBA case, Figure 36a, a valley extends across the 
space at = E-3.5. For the RR case the best second correction time de- 
creases to t^ = E-3. The values for the RR case are on the order of six 
times greater than the corresponding values of the RR + OBA case. 

At Uranus the best times for the second correction are again easy to 
identify for the two observation cases. For the RR + OBA case, Figure 37a, 
the best second correction time is about E-3.5. The RR case exhibits a nar- 
row horizontal valley at about t^ = 2 with a steep gradient above and 
below the valley. A vertical ridge extends across the space at about t^ = 
E-4.5. The values for the range-rate plus onboard angles case are 
roughly 12 times greater than the corresponding values for the RR case. 

5.4 Guidance Sensitivity Analysis 

The guidance survey presented in Section 5.2 provides important base- 
line data for preliminary assessment of the approach guidance propulsion re- 
quirements of the Grand Tour. A sensitivity analysis was performed on the 
nominal encounter profile at Jupiter, Saturn and Uranus to assess the sensi- 
tivity of these results to changes in the assumed error sources. For this 
analysis , the error sources and ground rules used for the reference navigation 
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accuracy analysis of Section 4.2 were retained. Both range-rate and onboard 
angle measurements were included. Additionally, the following approach guid- 
ance correction schedule was adopted: 

Jupiter - one correction at E-9 
Saturn - one correction at E-4 
Uranus - one correction at E- 3 . 

Results of the guidance sensitivity analysis are summarized in Table 6. 

Standard deviations of range-rate and onboard angle observation er- 
rors were decreased to one-tenth of their nominal values (90 percent decrease) 
to determine the effect of improved observation accuracy. As can be seen from 
iable 6, the approach guidance correction requirements are insensitive to the 
observation accuracy. Post-correction target error, however, is moderately 
sensitive to changes in both RR and OBA accuracy. For the 90 percent decrease 
in RR error the target error decreased 33, 36 and 23 percent at Jupiter, 

Saturn and Uranus respectively. The 90 percent decrease in OBA error resulted 
in an 80, 87 and 86 percent decrease in target error at Jupiter, Saturn and 
Uranus , respectively . 

One-sigma guidance correction errors were increased an order of mag- 
nitude (900 percent) to measure the effect of degraded guidance execution ac- 
curacy. Since only one approach guidance correction was specified, the ap- 
proacn velocity correction is not a function of guidance execution accuracy. 
Targe i error is relatively insensitive to guidance execution error. For the 
order of magnitude increase in the guidance magnitude errors , the resulting 
target error increases were 1.14, 2.26 and 7.90 percent at Jupiter, Saturn 
and Uranus, respectively. The corresponding target error increases due to 
increased orientation error were 1.72, 2.73 and 8.48 percent. 



Table 6 Approach Guidance Sensitivities 




Nominal 

Value 

Change in 
Parameter 

JUPITER, 

h = 

E-9 

SATURN, 

h = 1 

E-4 • 

URANUS, t 

1 = E 

- 3 

Error Source 

Parameter 

A V 

rms 

6T 

T 

rms 

AV 

rms 

<5T ' 
rms 

AV 

rms 





% 

m/ sec 

0 , 

0 

km 

Q, 

Q 

m/ sec 

% 

km 

% 

m/sec 

a 

0 


% 

Observations 
















Range-Rate 

a • 

P 

.5 mm/ sec 

-90 

0.00 

0.00 

-32 

-33 

0.01 

0.00 

-36 

-36 

0.00 

0.00 

-22 

-23 

Onboard Angles 

0 5 °o 

a g 

j 

10 arc-sec 

-90 

0.00 

0.00 

-78 

-80 

0.01 

0.00 

-86 

-87 

0.00 

0.00 

-83 

-86 

Guidance Correction 
















Magnitude 

a 

K 

.001 

900 

0 

0 

1.12 

1.14 

0 

0 

2.24 

2.26 

0 

0 

7.60 

7 . 90 

Orientation 

a 

Y 

.001 

900 

0 

0 

1.68 

1.72 

0 

0 

2.70 

2.73 

0 

0 

1 8.19 

8.48 

Initial State Error 

a s 

O 

* 

100 

4.88 

214 

0.10 

0.10 

15.19 

200 

0.15 

0.15 

19.1 

i 

200 

0.51 

9.52 


*See Table 1 
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The diagonal elements (variances) of the initial state error covariance 
matrix, S , were increased by a factor of four (100 percent increase in one- 
sigma values) to determine the effect of increased initial state errors. As 
can be seen from Table 6, the AV values increase by 214, 200 and 200 percent 
at Jupiter, Saturn and Uranus, respectively. Thus, is almost directly 

proportional to the standard deviation of the initial state errors. The post- 
correction target accuracy is relatively insensitive to the covariance of the 
initial state errors. 

5.5 Monte Carlo Simulation 

It is well known that the linear sequential estimation scheme used 
here yields poor results, and may diverge, if deviations from the nominal tra- 
jectory become. too large. Furthermore, the scheme can provide a satisfactory 
state estimate while yielding an unrealistic error covariance matrix for the 
estimate. The guidance algorithm is also based on linear theory and is in- 
valid for large state deviations. Therefore, it is necessary to test the va- 
lidity of guidance statistics predicted by the linear navigation and guidance 
theory. 

Initial State Errors 

A Monte Carlo simulation of the navigation and guidance process during 
Saturn approach was conducted for the purpose of checking the adequacy of the 
linear theory. The Saturn approach was selected because it is generally con- 
sidered to be the most critical encounter from a guidance standpoint. A one- 
correction guidance logic was adopted and a simulation was carried out for 
two different correction times: t^ = E-10 and t = E-4. For this analysis, 

the error sources and ground rules used for the guidance survey were retained 
except that the ’true" Initial state vector was simulated by adding a random 



error vector to the nominal state vector. Components of the error vector 
were assumed to be uncorrelated and each was generated by sampling from its 
respective normal distribution (see Section 2.3). The distribution for each 
component is defined by a zero mean and a variance equal to the diagonal ele- 
ments of the covariance matrix S . Since S = P , the variances are those 

o o o 

values listed in Table 2. Values of the mean and standard deviation of each 
component are listed in Table 7 for the 35 simulation runs. 

Monte Carlo Results 

The results of the simulations are presented in Figure 38. Actual 
values of the target parameters B*R and B-T are plotted with respect to the 
nominal target point. The standard probability ellipses representing the 
bivariate distribution of the predicted target error are shown for compari- 
son. Pertinent statistics of the simulated data are contained in Table 8. 

The table contains values for the mean and standard deviation of the target 
errors, the magnitude of the velocity corrections, AV ? , and the magnitude 

O O ^2 

of the error vector, 1 6 T T | = [6B*R + SB-T J . 

Even for the relatively small sample size (35 runs), the correlation 
between the simulation statistics and the predicted statistics is reasonably 
good for the t^ = E-10 case and very good for the t^ = E-4 case. For the 
E-10 case, Figure 38a, the mean of the simulation errors is 2.3 km for 6B«R 
and is -14.0 km for 6B-T compared to predicted values of zero for both. 
Standard deviations of 6B*R and 6B*T are 185 and 82 km, respectively, com- 
pared with predicted values of 201 and 92 km. Theory predicts a AV 1 RMS 
value one percent greater than the value obtained by simulation (3.03 m/sec 
compared to 3.00 m/sec). The actual RMS value for the magnitude of the tar- 
get error, | 6T T | is 203 km compared with a predicted value of 221 km. 
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Component 

x (km) 
y (km) 
z (km) 
x (km/day) 
y (km/day) 
z (km/day) 
b^ (km) 
b (km) 

y 

b (km) 
z 


Table 7 Simulated Initial State Error Statistics 


Mean 


One -Sigma 


-122 

-196 

204 

3.06 

3.33 

-1.18 

-350 

108 


1432 

1010 

1473 

16.2 

13.5 

16.7 

1466 

1118 


93 


1288 
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400 1 
AB-R 

(38b) E-4 




Figure 38 Monte Carlo Simulation Results 
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For the E-4 case. Figure 38b, the mean of the simulation errors is 
2.7 km for 6B-R and is 1.2 km for 6B*T compared to predicted values of zero 
for both. Standard deviations of 6B*R and 6B*T are 87 and 46 km, respec- 
tively, compared with predicted values of 86 and 49 km. The AV f RMS value 
of the simulation is 7.7 m/sec compared with 7.6 m/sec predicted by theory 
(1.5 -6 greater). The actual RMS value for the magnitude of the target error 
is 109 km compared with a predicted value of 99 km. 

Correlation between theory and simulation should improve as the sam- 
ple size is increased. However, based on the analysis discussed above, it 
can be concluded that the navigation and guidance algorithm formulated in 
Chapter 2 provides an adequate model for preliminary mission analysis of the 
Grand Tour. For the example selected for simulation, the theory slightly 
over-predicted the simulation RMS values of |6T T | and AV T . For the E-10 
case, RMS values of AV T and | 6T 1 | are 6 and 17 percent larger, respectively, 
than the corresponding mean values. For the E-4 case the RMS values of AV* 
and |6i f | are 7 and 13 percent larger than the corresponding mean values. 
Thus, RMS values of correction magnitude and target error magnitude appear 
to be good first approximations to the mean of their non-Gaussian distribu- 


tions . 



CHAPTER 6 


CONCLUSIONS AND RECOMMENDATIONS 

Ihe investigation reported here was concerned with the simulation 
and analysis of the navigation and guidance process for the Jupiter, Saturn 
and Uranus encounter phases of the 1977 Grand Tour Interior Mission. 

Principal Conclusions 

From an examination of the theoretical development of the navigation 
and guidance algorithms in Chapter 2, it can be concluded that the error in 
the position of the outer planets can be conveniently included as a parame- 
ter to be estimated if the error is assumed constant over the encounter 
period. 

Based on a study of the navigation data presented in Chapter 4 the 
following conclusions can be drawn: 

1. Range-rate data provide adequate in-plane navigation accuracies 
but poor out-of -plane navigation accuracies during the approach 
phase at each planet . 

2. Although onboard angle data provide relatively poor navigation 
accuracies in terms of the planetocentric state of the space- 
craft 5 the data yield good estimates for the position components 
of the target vector, i.e., B-R and B-T. 

3. Combining onboard angle data with range-rate data leads to a sig- 
nificant improvement in overall navigation accuracy at each 
planet . 
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4. Range data, when combined with range-rate data, provide only 
slight improvement over the range-rate only data at each planet. 

5. Onboard angle data alone does not significantly reduce planeto- 
centric spacecraft position uncertainty at Jupiter until the 
spacecraft reaches periapsis. At Saturn onboard angle data yield 
better planetocentric position accuracy than do the range-rate 
data up to about E-16 days. The same is true at Uranus where the 
crossover point is about E-10 days. Near periapsis the range- 
rate data is superior for all three planets . 

6. Range-rate data provide better planetocentric velocity accuracy 
during planet approach than do the onboard angle data. Combin- 
ing onboard angle data with range-rate data does not result in a 
significant decrease in velocity uncertainty over the range-rate 
only case except at Uranus. 

7. Range-rate observations alone provide a significant improvement 
in Jupiter and Saturn position uncertainty, but only a slight im- 
provement in Uranus position uncertainty. Onboard angle data 
alone provide virtually no planet position information. 

8. Combined range-rate and onboard angle data provide a signifi- 
cantly smaller planet position uncertainty than the range-rate 
data alone. 

9. Range-rate data provide practically no B-R and B-T information 
until the spacecraft is well within the planet T s sphere of in- 
fluence. As the spacecraft nears encounter the B-T uncertainty 
decreases rapidly while the B*R uncertainty decreases at a much 


slower rate. 
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10. Onboard angle data provide good B*R and B*T accuracy but poor t 
accuracy . 

11. Adding star-planet onboard angle (an out-of-plane angle) observa- 
tions to the range-rate observations decreases B-R uncertainty 
significantly. A combination of all three observation types 
(range-rate plus both onboard angles) yields good target accuracy 
in all three target components . 

From a study of the guidance analysis results presented in Chapter 5 
the following conclusions are reached: 

1. A survey of approach AV requirements for a two-correction 

rms 

guidance logic indicates that: (a) the sum of the approach and 

departure AV values is sensitive to the time of the second 
rms 

correction, but relatively insensitive to the time of the first 
correction; (b) the sensitivity to the second correction time is 
greater for the range-rate case than for the range-rate plus on- 
board angle case; (c) the optimum correction times for the second 
approach corrections is later for the range-rate plus onboard 
angle case than for the range-rate case, and (d) the optimum cor- 
rection times for the second approach correction are near the 
following survey values: 



RR 

RR + OBA 

Jupiter 

E-9 

E-10 

Saturn 

E-3 

E-3. 5 

Uranus 

E-2 

E-3. 5 


2. A sensitivity analysis of the navigation and guidance process 

demonstrated that: (a) post-correction target error is relatively 
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insensitive to both initial state error and guidance execution 
error; (b) post-correction target error is moderately sensitive 
to both range-rate and onboard angle observation accuracy, being 
more sensitive to onboard angle accuracy than to range-rate ac- 
curacy; and (c) guidance correction requirements are insensitive 
to observation accuracy but are almost directly proportional. to 
the standard deviation of the initial state errors . 

3. Results of a Monte Carlo simulation of the navigation and guid- 
ance process verify that the linear estimation and guidance 
theory provides an adequate model for preliminary mission analy- 
sis of the Grand Tour. 

Recommendations for Future Study 

The research reported here uncovered several areas which should be 
considered for future study. The following specific studies are recommended: 

1. An extension of the analysis reported here to include other ob- 
servation types (e.g., moon-star and moon-planet angles) and 
other missions (e.g., the three planet Grand Tour missions). 
Accordingly, the computer programs STEP V and STEP VI should be 
modified to include the capability to simulate and process addi- 
tional observation types. 

2. A study of the effect of other error sources (e.g., unmodeled 
accelerations, planet mass uncertainty, tracking station loca- 
tion errors, and planet center finding bias) on navigation accu- 
racy and guidance requirements . Computer programs STEP V and 
STEP VI could be modified to provide options for estimating (or 
including but not estimating) additional error parameters. 
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3. A study of the effect of secondary perturbations -on the motion 
of the spacecraft during encounter. A more sophisticated dy- 
namical model which includes planetary oblateness and the per- 
turbation of additional bodies should be considered for incor- 
poration into the STEP V and STEP VI computer programs. 

4. A study of alternate guidance parameters for planetary flyby 
missions should be conducted to determine the best target vec- 
tor for the Grand Tour missions . 

5. The interactive graphics feature of computer program STEP VI 
should be further developed to exploit the full capabilities of 
the CDC 6600/252 interactive system. For example, additional 
problem status displays could be added to automatically display 
pertinent data as problem execution proceeds. Also, use of the 
light pen feature would greatly improve the manipulation of in- 
put data. 

6. An analysis of information lag times due to light travel time 
should be made to determine the effect on navigation and guidance 
requirements . 
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APPENDIX A 

NOMINAL GRAND TOUR MISSION DATA 


Mission/traj ectory data for the 1977 Grand Tour Interior mission 
are presented in the following tables. Trajectory parameters which define 
the four heliocentric legs of the mission are shown in Table 9. The classi- 
cal elements which define the orientation of the transfer conics (i, 8, w) 
are defined with respect to the equinox and ecliptic of date. Parameters 
which define the encounter trajectory are presented in Table 10. The clas- 
sical elements which define the orientation of the encounter conics are de- 
fined with respect to a planetocentric equatorial coordinate system. 

The nominal mission data were generated using a computer program 
developed by M. C. Poteet of Convair Aerospace, Fort Worth, Texas. The pro- 
gram generates an approximate solution to multi-planet flyby missions by 
solving for the heliocentric transfer conics which match the hyperbolic ex- 
cess speeds at arrival and departure. The heliocentric solution defines 
heliocentric velocity requirements , and thus , the hyperbolic asymptotes , at 
arrival and departure. The asymptotes and the excess speed completely de- 
fine the planetocentric flyby conic. 
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Table 9 Heliocentric Trajectory Parameters 


LEG 

JULIAN DATE 244- 

DURATION 

a 

e 

i 

n 

w 

TRANSFER 


DEPART 

ARRIVE 

DAYS 

10 6 km 


deg ■ 

deg 

deg 

ANGLi. 

deg 

Earth-Jupiter 

3390.5 

3895.8 

505.3 

951.9 

.842 

.867 

-19.11 

2.38 

145 

Jupiter-Saturn 

3895.8 

4472.2 

576.4 

-459.2 

2.690 

2.449 

111.18 

- .206 

43 

Saturn-Uranus 

4472.2 

5701.8 

1229.6 

-258.6 

6.470 

2.379 

70.74 

113.34 

76 

Uranus -Neptune 

5701.8 

6700.0 

998.2 

-208.6 

12.151 

2.365 

-112.28 

“34.98 

62 



Table 10 Planetocentric Flyby Parameters 





ENCOUNTER CONDITIONS 

EXCESS 





TAR/eFT 


| ENCOUNTER DATE 

ALTITUDE 

VELOCITY 

VELOCITY 

<$> 

a 

10^ km 

e 

i 

PARAMETER, b 


JD 244- 

CIVIL 

103 km 

RADII 

km/sec 

km/sec 

deg 


deg 

Jupiter 

3895.5 

22 Jan 79 

194.5 

2.73 

33.23 

12.3 

98.8 

-838.4 

1.32 

5.6 

718.7 

Saturn 

4472.2 

20 Aug 80 

2.4 

0.04 

38.62 

16.9 

85.5 

-132.7 

T; 47 

29.2 

143.6 

Uranus 

5701.8 

2 Jan 84 

13.1 

0.56 

27.85 

21.4 

29.7 

- 12.6 

3.90 

100.1 

47,6 

Neptune 

6700.0 

26 Sep 86 

OPEN 

23.9 



op: 

EN 


— — - — 1 
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THE DYNAMICAL EQUATIONS IN RECTANGULAR CARTESIAN COORDINATES 


The differential equation of state, Eq. (2.9), can be written in a 
rectangular cartesian coordinate system as 


where 



F„ = 


F„ = 


F„ = 


F. = - 


yx 

- y s 

1 

X 

+ x t 


3 

r 

- 

CO 

04 

X 

py 


y 

+ Y. 

L 

V 

3 

r 

" V 


CO 

04 

X 

pz 


z 

+ z t 

V 

CO 

1 fc 
1 
1 

■ P s 


CO 

04 

X 


(B.l) 
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where X t *, Y t *, and Z ± * are the components of the nominal heliocentric posi- 
tion vectors obtained from the planetary ephemeris subroutine and b , b and 

x y 

b z are components of bias in the planet position vector. The bias in the 
planetary position vector represents an error in the ephemeris of the target 
planet which is assumed to be constant over the time period of interest. 

The state transition matrix is evaluated by numerical integration of 
the matrix differential equation 


3>(t,t Q ) = A(t ) 4>(t,t ) 


(B .3) 


with initial conditons $>(t ,t ) = I 

o o ’ 

is a 9x9 matrix of partial derivatives evaluated on the 
nominal trajectory. For convenience the A-matrix can be partitioned into 
nine submatrices as follows : 


where A(t) 


df 

JX_ 



(B.4) 
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The symmetric submatrices A^ and are defined below. 





121 


5 

J 47 3b 


= U 


x 


3X 

5 


3X, 


R 


-1 + 1 

R 3 R 3 
t 


- 5 

i 48 3b 


y 


3XY 

5 


R 


3X t 2 t 
R 5 


- !!i 

*49 3b 


= 4 


3XZ 

„5 


3X t Z t 


3F, 


57 3b 


48 


^ F c 

1 58 " 3b“ 


= 4 


y 


3Y 

R 5 


3Y, 


i + i_ 

R 3 R. 3 


d Jj 

*59 " W 


= 4 


3YZ 

5 


R 


3Y t Z t 
R. 5 


. .5 

67 3b 


= a 


x 


49 


3F, 


68 3b 


y 


59 


, =5 

69 db 


= U 


3Z" 

5 


3Z. 


R 


i + 1 

R 3 R 3 
t 


(B.6) 
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the observation-state relationships in rectangular 


CARTESIAN COORDINATES' 


functional relationships between the observations and the compon- 
ents of the state vector can be written 
system as follows : 


m a rectangular cartesian coordinate 


where 


P - I(x - x ) 2 + (y - y ) 2 + (z _ z ) 2 i 

S S 1 

P = [(X - X )(x - x ) + (Y - Y )( 

b s s 

“ = cos ’L<>-.X t + yY t + Z z t t r 2 )/(rR)J 


s Y - V < 4 


- 1 , 


8 cos {-(x s + ys + ZS )/rJ 

x y 2 


(C.l) 

Z s )(Z-Z s )]/p (c.2) 

(C.3) 

(C . 4) 


x t * + b^ + X 

(C . 5 ) 

Y * + b + Y 

t y 

(C.6) 

Z * + b + z 
t z 

(C . 7 ) 

[X 2 + Y 2 + Z 2 ]^ 

(C.8) 

[x 2 + y 2 + z 2 J 2 

(C . 9 ) 


and S x’ S y ’ and s z ars the direction cosines of the navigation star direction. 
The linearized observation-state relationship 


y(t) = H(t ) x ( t ) 


(C.10) 


m component form becomes 
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"6P 


h ll h 12 

h 19 


Sx 

6p 

__ 

h 21 



6y 

6a 


h 31 



6z 

JL 


_ h 41 

49 


6u 


<$v 

6w 


6b 

x 

6b 

y 

6b 

L 3 

where [ J* indicates that the elements of the H(t) matrix are evaluated on 
the reference trajectory, X*(t). 

r\j 

The elements of the H-matrix are as follows : 


h n ■ fe = - V /p 


12 


5p 

3y 


= (Y - 


v s )/» 


h 13 " 3z ' (Z 


2 S )/P 


, 3p 

h 14 = H = 0 


15 


5p 

3v 


16 


Sp_ 

Sw 


= 0 
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h - is: - v. 

17 3b 11 

x 

h - iL - h 

18 3b 12 

y 

h 19 3b h 13 

z 

h 2i = 3 - V - (x - v (p/p > 1/p 

h 22 = F = I(X - V - (Y - Y s )(p/p)]/p 

h 23 = || [4 - Z s ) - (Z - Z £ )(p/p)]/p 

h = 3fi. - 

h 24 3u h ll 

h = li = h 

25 3v 12 


h 


26 


3p_ _ h 

3w 13 


h 


27 


3p. _ h 

3b 21 


h 


28 


ap. _ h 

3b " 22 

y 


h 


29 


3p _ h 

3b 23 


h 


31 


da 

dx 



r x 

xl 

cos a 

2 _2 

— 

Lr 

R -J 


(x+X) 

rR 


1 


sin a 
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46 
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3w 


h 36 _ n 

h 47 3b “ ° 

x 


h - 36. - n 

h 48 3b ° 

X 


h 49 If = 0 
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APPENDIX D 


THE DERIVATION AND APPLICATION OF THE SCHUR IDENTITY 


The Schur identity (also called the "inside out rule") can be stated 
as follows : 

The matrix equation 

E = (A + BC)" 1 (D.l) 

is equivalent to 

E = A -1 - A -1 Bfl + CA^BJ 1 CA _1 (D.2) 


where A and E are nonsingular nxn matrices and where B and C are any matrices 
with dimension nxp and pxn respectively. 

Derivation j 

Given 

E = (A + BC) _1 (D. 3) 

then 

E _1 = A + BC . (D.4) 

Post-multiplying by E gives 

I = AE + BCE (D.5) 


and pre-multiplying by A 


yields 


A 1 = E + A" 1 BCE 
E = A"" 1 - A -1 BCE 


(D.6) 
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Now to eliminate E from the right hand side. First pre-multiply Eq. 
(D . 7 ) by C 

CE = CA" 1 - CA _1 B CE 
CE + CA"' 1 BCE = CA - ' 1 

(D . 7 ) 

[I + CA -1 B]CE = CA -1 • 

CE = II + CA -1 B.]~ 1 CA -1 

Then substitute Eq. ( D - 7 ) into Eq. (D.6) to obtain the desired result: 

E = A -1 - A -1 BlI + CA -1 BJ _1 CA -1 . (D.8) 


Application 


Consider the nxn. covariance matrix 


P = (.P 1 + H I R H ) 

, k k k k k 


-1 


(D.9) 


which is of the form 


E = (A + BC) 


-1 


_ -1 T -1 

where A = P^ ' , B = , and C = R ' H^. Applying the Schur identity (Eq. 

D.8) gives 


p k = f k - ^ * R 'Vk H kJ~ R_1 H k f k • 


(D.10) 


Eq. (D.10) can be rewritten to obtain the desired result, i.e.. 


p k = \ - K p k H k T + R)_1 H k ? k • 


(D.ll) 
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THE DERIVATION AND APPLICATION OF THE ELLIPSOID OF PROBABILITY 


In the analysis of navigation and guidance errors, it is sometimes 
convenient to display in geometric form the statistical data contained in a 
covariance matrix by generating surfaces of constant probability - 

Consider a multivariate Gaussian distribution of the random n-vector 
z with zero mean and covariance 


A = E[z z T ] . 

The distribution is characterized by the quadratic form 


(E . 1 ) 


a = 


T A -1 
z A 


z - 


A + zz 


- 1 


(E.2) 


The scalar q can be shown to have a chi-square distribution with zero mean, 

(19) 

unit variance and n-degrees of freedom . Setting q = q^ in Eq. (E.2) de- 
fines an "ellipsoid of probability." The probability that z falls within the 
ellipsoid corresponding to q^ is the probability that q < q^ and is given by 

i f qi t 2 i i> 

P (q ,n) = —pr q exp(-q/2)dq (E.3) 

r 1 2 n/ i (n/2 ) Jo 

The lengths of the principal axes of the error ellipsoid are given by 


a. 

i 



A. . 

i 5 


i = 1, . . . ,n 


(E.4) 
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where A^ is the i— eigenvalue of the covariance matrix A. The orientations 
of the principal axes are given by the unit vectors 

X. 

e i = TxTT ; i = 1, . . . ,n (E. 5 ) 


where is the i— eigenvector of the covariance matrix A. 


The Bivariate Case 


Consider the 2x2 covariance matrix 


A = 


the eigenvalues of A are formed from 


a ll 

a l2 

a J 

a 

21 

22 


(E.6) 


A - AI = 0 


(E.7) 


In expanded form, Eq. (E.7) becomes 


a . . - A a . 0 
11 12 


a 21 a 22 X 


X 2 - X(a 1;L + a 22 ) + a n a 22 a 12 a 2 i 


Now let 


2A ( a 11 + a 22 ) 1 + a 22' )2 4< ' a il a 22 “ a l2 a 21^ * 


= rr2 


a H = a 


a 22 °y 


S 12 a 21 


a = pa a 
xy x y 
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Then 5 


2A = ( a 2 + a 2 ) ± {(a 2 + a 2 ) 2 - 4(a 2 c 2 - p 2 o 2 a 2 )]^ 

x y xy xy xy 

2A = (a 2 + a 2 ) ±. [a 4 + 2a 2 a 2 + a 4 - 4a 2 a 2 + 4p 2 a 2 a 2 ]^ 

x y x xy y xy xy 

2A = (a 2 + a 2 ) ± -[(a 2 - a 2 ) 2 + 4a 2 ] 2 

x y x y xy 


which yields the eigenvalues 


(a 2 + a 2 ) 
x 


1,2 


± [^(a 2 - a 2 ) 2 + a 2 


xy 


(E.8) 


The two eigenvectors = 

X 1 

_ y l_ 

and X^ = 

X 2 

_ y 2 

or 




[A - A i IJX jL 

= o ; 

i = 


are defined by AX. = A.X.: 

J i li 


(E . 9 ) 


Hence from Eq. (E.9) 


a . „ - A . a„ „ i 


F— 

x. 

11 l 12 1 


1 

_ a 21 a 22 " A ij 


_ y i_ 


U 11 - V X i + a 12 y i = ° 


a 21 X i + (a 22 ‘ = 0 


Setting x^ = 1 gives 


(a U - V + a !2 y i = 0 


A. - a . . 

l 11 


12 


i = 1,2 


(E.10) 


The eigenvectors are then 
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i = 1,2 . (E.ll) 


The orientation of the error ellipse is specified by the angle <j> be- 
tween the major axis of the ellipse and the x-axis. The angle is given by 


cos <j> = e^ • i 
sin c/> = e 2 - j 

where 



(E.12) 


This can be reduced to 



(E . 13 ) 


The semimajor and semiminor axes of the error ellipse are given by 


a 


b 


I* A 




For the bivariate case Eq. (E.3) reduces to P (q^,2) = 1 - exp(-q /2). 
The probabilities corresponding to /q^ - 1,2, and 3 are 0.40, 0.86, and 0.99 
respectively. Values of q^ corresponding to other given values of P_^ can be 
computed from q^ = 2£n(P^ - 1). 


Target Dispersion Ellipse 

Consider the error covariance matrix associated with the target 
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coordinates B*R and B*T 


A = 


B-R 


B*R ? B*T 


B«R,B*T B*T 


Figure 39 shows the standard probability ellipse (q = 1, = 0.40) which 

represents A. The ellipse is defined by 


a = 


I a 2 + 

‘ B-R B-T 


2 - + ^B-R + + °B.R,B.F 


v I ^ 
" 5 ' 


l0 B • R + a B-T n , 2 ^ 2 n 2 ^ 2 i"2 

2 ^ a B*R + °B-T + a B*R,B’T 


id h 


= tan 


a 2 - a 2 
-1 I B-R 


B-R, B-T 


= 90 + 4> 


-90 < <j> < 90 
0 < ip < 180 
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DEFINITION OF THE TARGETING COORDINATE SYSTEM 


The targeting vector used in this study is defined by 


B * R 1 


T 


B-T 


t 

e 


(F . 1 ) 


where t^ is the time of encounter (periapsis) and where B-R and B-T are the 
components of the vector B in the target (R - T) plane. The vector B, de- 
fined as the normal vector directed from the center of the target planet to 
the approach asymptote of the encounter hyperbola , is illustrated in Fig. 
40 . 


R - S - T Coordinate System 

Let S be the unit vector in the direction of the approach asymptote. 
The unit vector T is orthogonal to S and parallel to the reference plane 
(e.g., ecliptic plane), i.e., 


T = 


(S MP ref )/| 


S X P 


ref 1 


(F . 2 ) 


where P ^ is the pole of the reference plane. The unit vector R completes 
ref 

a right-handed orthogonal system, i.e.. 


R = S x T . (F . 3 ) 

The B- vector targeting system involves the intersection of three planes, i.e., 
the orbit plane, the reference plane, and the target plane. The following 
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Figure 40 Targeting Geometry 
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notes may be helpful in visualizing the system: 

(1) T lies along the intersection of the reference plane and the 
target plane. 

(2) B lies along the intersection of the orbit plane and the target 
plane . 

(3) The orbit plane and the target plane are orthogonal. 
Computational Algorithm 

Given the state of a spacecraft in hyperbolic orbit relative to the 
planet (i.e., r and v) and the gravitational parameter p, the vector B is 
determined as follows : 


1) 

r = 

(r-r)' 2 

distance 

2) 

v = 

(v-v) 3 * 

speed 

3) 

C 3 

2 

= v - 2y/r 

energy constant (twice the 

4) 

V 

00 


energy per unit mass) 
hyperbolic excess speed 

5) 

w = 

(r x v)/ |r x v | 

unit vector normal to orbit 

6) 

h = 

v- (W x r ) 

plane 

magnitude of angular momentum 

7) 

a = 

U /C 3 

semimajor axis 

8) 

b = 


semiminor axis 

9) 

P = 

(v 2 /y-l/r)(l/e)r - (r-v/y )(l/e)v 

unit vector in the direction 

10) 

Q = 

W x P 

of periapsis 

unit vector normal to W and P 

11) 

B = 

aCe-l/e )P - aA-l/e* Q 

target vector 

12) 

H = 

_ ^ 

cos [ae/(a+r )] 

hyperbolic auxiliary angle 

13) 

t 

go 

3 3^ 

= (H /p) 2 {e tan H - Logftan(rr/4 + 

H/2 ) ] } 

time to encounter (periapsis) 

14) 

t = 
& 

= t + t 

time of encounter 
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Additional parameters of the hyperbolic orbit may be computed as follows: 


15) 

r 

P 

= a(e-l )P 

periapsis radius vector 

16) 

= 

-1 

it -2 tan (b/a) 

turn angle 

17) 

f = 

-2 tan h/(e+l)/(e-l) tan(H/2)] 

true anomoly 

18) 

S = 

(W x B )/ | W x B| 

approach asymptote 

19) 


= cos ip S - sin ip (B/|b|) 

departure asymptote direction 

20) 

i = 

cos ^(W ) 
z 

inclination 

21) 

a = 

tan _1 (W /-W ) 
x y 

ascending node 
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DERIVATION OF THE TARGET MAPPING MATRIX C 


Consider the departure target vector expressed as a function of 

°°D 

the approach target vector T by a Taylor expansion about the nominal path. 


i . e . , 


CO CO 

D D 


3V 

oo 



6T + 


(G.l) 


Neglecting terms higher than first order yields 


where C 


3V 

OO 

D 

3T 


6V = C <5T 

OO 

D 


(G.2) 


In what follows an analytical expression for the linear mapping matrix C 

will be derived for a two-body dynamical model. 

Consider the geometric relationship between the departure target 

vector , the approach asymptote S 5 and the miss vector B shown in the 
°°D 

sketch below. 
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The departure target vector can be written as 


V = V [cos tp S - sin ip — ] 

oo oo T h J 

D D D 


where ip is the turn angle, b = |B|, and = |v |. Resolving B into the 

"D 

targeting components B*R and B-T yields 


= V^Icos ip S - sin ^{(B-R)R + (B*T)T}] . 
D 


Taking — gives 


3V 

OO 

D 

3B-R 

V 

JT< b 


simp 


b 


simp 


b 

3V 

oo 

D = 
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/V \ 

(-)» 


simp 


b 


simp 


b 

3V 

OO 

D 

at 

e 

0 

2 

+ B-T , 

then 


{sin \p + (B-R) 


B .R } R - { f|±g* (B . T) 


B • T}T] 


Scosip - r . 3sinip \ 

dB — s - {sin * + WIT (B ' R) 


}R - ll§if (B'T) 


(*>— \ 
\3B-R ) 

l—\ 

\dB-R) 


B • T}T] 


3b _ B-R 
3B-R b 
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and 


3b 

3B * T 


B-T 


The turn angle is defined by 


cos b - 


2 v 4 \ 

- 1 


b V 


/ bV \ 


y 2 


I 


+ l 


\ / 


thus , 


dCOSlp 

3b 


4b V 


2 4 \ 

b V ' 

■2 1 t 1 


and 


Define 


Bsinij; / 9cosi|A , , v 

Sb - = (lb— 


and 


Then 



A - 

4 

4b V 

00 


d l " 

1 b 2 V 4 


y 2 | 

00 

l “ 2 


d 2 = ~ d i cot 

3cos ip 

3cosi(j 

db 

3B-R 

3b 

dB-R 

3cos^ 

3cos<jj 

db 

3B-T 

3b 

dB-T 


+ 1 


= d. 


B-R 
1 ( b 


= d. 
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Likewise 
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Now Eqs. (G 
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and defining 
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leads to 
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